Computation of gravitational waves from inspiraling binary neutron stars in 
quasiequilibrium circular orbits : Formulation and calibration 



1 Masaru Shibata and 2 Koji Uryu 
1 Graduate School of Arts and Sciences, University of Tokyo, 
Komaba, Meguro, Tokyo 153-8902, Japan 
2 Department of Physics, University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI53201, USA 



o 
o 

(N 
Or 



> 

O 

o 



o 

i 

5h . 

bo 



Gravitational waves from binary neutron stars in quasiequilibrium circular orbits are computed 
using an approximate method which we propose in this paper. In the first step of this method, 
we prepare general relativistic irrotational binary neutron stars in a quasiequilibrium circular or- 
bit, neglecting gravitational waves. We adopt the so-called conformal flatness approximation for 
a three-metric to obtain the quasiequilibrium states in this paper. In the second step, we com- 
pute gravitational waves, solving linear perturbation equations in the background spacetime of the 
quasiequilibrium states. Comparing numerical results with post Newtonian waveforms and luminos- 
ity of gravitational waves from two point masses in circular orbits, we demonstrate that this method 
can produce accurate waveforms and luminosity of gravitational waves. It is shown that the effects 
of tidal deformation of neutron stars and strong general relativistic gravity modify the post New- 
tonian results for compact binary neutron stars in close orbits. We indicate that the magnitude of 
a systematic error in quasiequilibrium states associated with the conformal flatness approximation 
is fairly large for close and compact binary neutron stars. Several formulations for improving the 
accuracy of quasiequilibrium states are proposed. 
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I. INTRODUCTION 

The last stage of inspiraling binary neutron stars to- 
ward merger, which emits gravitational waves of fre- 
quency between ~ 10 and ~ 1000Hz, is one of the most 
promising sources of kilometer-size interferometric grav- 
itational wave detectors such as LIGO Q]. Detection of 
gravitational waves from the inspiraling binaries will be 
achieved using a matched filtering technique in the data 
analysis, for which it is necessary to prepare theoretical 
templates of gravitational waves. This fact has urged 
the community of general relativistic astrophysics to de- 
rive highly accurate waveforms and luminosity of gravi- 
tational waves from compact binaries. 

For an early inspiraling stage in which the orbital sep- 
aration r a is > 4i? where R denotes neutron star radius 
and in which the orbital velocity v is much smaller than 
the speed of light c, tidal effects from companion stars 
and general relativistic effects between two stars are weak 
enough to neglect the finite-size effect of neutron stars as 
well as to allow us to adopt a post Newtonian approxi- 
mation. For this reason, post Newtonian studies jointly 
using point particle approximations for compact objects 
have been carried out by several groups, producing a wide 
variety of successful results (e.g., [0-0). However, for 
closer orbits such as for r < 4i? and v > c/3, the tidal 
effect is likely to become important, resulting in deforma- 
tion of neutron stars and in the modification of the ampli- 
tude and luminosity of gravitational waves. Furthermore, 
general relativistic effects between two stars are so signif- 



icant that convergence of post Newtonian expansion be- 
comes very slow 0. These facts imply that, for preparing 
theoretical templates for close orbits, fully general rela- 
tivistic and hydrodynamic treatments for the computa- 
tion of binary orbits and gravitational waves emission are 
necessary. 

Using the quadrupole formula of gravitational wave lu- 
minosity dE/dt and the Newtonian formula for the bind- 
ing energy between two point masses, E p , the ratio of co- 
alescence timescale due to emission of gravitational waves 
E p / (AdE/dt) |)| to the orbital period for binaries of equal 
mass in circular orbits is approximately written as 
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where M t is the total mass and G gravitational constant. 
The effects of general relativity and tidal deformation can 
shorten the coalescence timescale by a factor of several 
(see Sec. V), but for most of close orbits, the emission 
timescale is still longer than the orbital period. This 
implies that binary orbits may be approximated by a 
quasiequilibrium circular orbit, which we here define as 
the orbit for which the coalescence timescale is longer 
than the orbital period. 

Several approximate methods with regard to the com- 
putation of the quasiequilibrium states and associated 
gravitational waves have been recently presented by sev- 
eral groups ]lO|dTz|. All these methods require one to 
solve the Einstein equation by direct time integration and 
hence require to perform a large-scale numerical simula- 
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tion. In 00, new formalisms have been proposed to 
compute the late inspiraling stage of binary black holes, 
for which it is likely to be necessary to perform numerical 
time integration for obtaining a realistic quasiequilibrium 
sequence. On the other hand, the purpose for the authors 
of |l2"f is to compute gravitational waves from a fixed 
background spacetime of a computable quasiequilibrium 
state such as that of binary neutron stars. However, so 
far, the two formalisms have not been applied yet pO| , pl| , 
and the other method has not succeeded in an accurate 
computation of gravitational waveforms because of re- 
stricted computational resources iQ. To adopt these 
methods in accurately computing gravitational waves, it 
is necessary to develop robust computational techniques 
as well as to prepare sufficient computational resources 
for a large-scale simulation. 

The purpose of this paper is to compute gravita- 
tional waves from binary neutron stars in quasiequilib- 
rium states. A quasiequilibrium sequence of binary neu- 
tron stars can be constructed characterizing the sequence 
in terms of conserved quantities such as baryon rest mass 
and vorticity. In addition, some approximate formula- 
tions and numerical techniques have been already devel- 
oped for the computation of such quasiequilibrium solu- 
tions |l^-[T7[|. These facts imply that we may avoid per- 
forming direct time evolution of the Einstein and hydro- 
dynamic equations for obtaining binary neutron stars in 
quasiequilibrium. Only in computing gravitational waves 
do we need to integrate the Einstein equation using the 
quasiequilibrium solution as a source. From these rea- 
sons, we follow an idea of |12| ], but we propose a more 
systematic approximate formalism in which it is possible 
to compute waveforms and luminosity of gravitational 
waves with better accuracy using well-known computa- 
tional techniques and cheap computational costs. 

Our method is in a sense similar to the standard post 
Newtonian method for the computation of gravitational 
waves from binaries of two-point masses in circular or- 
bits &§]. Thus, before proceeding, let us briefly re- 
view an outline of the post Newtonian method. In the 
post Newtonian calculation, the procedure is divided into 
two steps: In the first step, the quasiequilibrium circular 
orbits of binaries are determined using post Newtonian 
equations of motion for two point masses, neglecting ra- 
diation reaction terms of gravitational waves. Neglect 
of the radiation reaction is justified for most of orbits 
for which the radiation reaction timescale is longer than 
the orbital period as shown in Eq. (1.1). After the bi- 
nary orbits are determined, gravitational waves are calcu- 
lated in a post-processing; one integrates the post New- 
tonian wave equations for gravitational waves, substitut- 
ing the matter field and associated gravitational field of 
quasiequilibrium states as the source terms. After the 
computation of the gravitational wave luminosity, one 
can compute the radiation reaction to a quasiequilibrium 
circular orbit to determine a new orbit. By repeating this 
procedure, one can determine an evolution of a binary or- 
bit due to radiation reaction of gravitational waves and 



associated gravitational wave train. 

As in the post Newtonian method, in our formalism, 
quasiequilibrium states are computed in the first proce- 
dure assuming that gravitational waves are absent. As 
a first step of the development of our new scheme, we 
adopt the so-called conformal flatness approximation for 
computation of the quasiequilibria in this paper. After 
computation of the quasiequilibrium states, we integrate 
the wave equation for gravitational waves (derived from 
the Einstein equation in Sec. IV), inputting the gravi- 
tational and matter fields of the quasiequilibrium states 
as the source terms. The difference between the post 
Newtonian method and ours is that we fully take into 
account general relativistic effects (under the adopted 
approximate formulation) and hydrodynamic, tidal de- 
formation effects. As is shown later, these two effects 
play important roles for compact binary neutron stars in 
close orbits. 

A word of caution is appropriate here: We choose 
the conformal flatness approximation for the quasiequi- 
librium solutions simply because of a pragmatic reason 
that we currently adopt this approximation in numerical 
computation. It would be possible to extend this work 
modifying the formalism for the gravitational field of the 
quasiequilibrium background solutions (see discussion in 
Sec. VI). The purpose in this paper is to illustrate the 
robustness of our new framework. 

The organization of this paper is as follows. In Sec. 
II, we describe the Einstein equation in the p resen ce of 
a helical (helicoidal) Killing vector [cf. Eq. (2.1)]. In 
deriving the equations, we do not consider any approxi- 
mation and assumption except for the helical symmetry. 
We will clarify the structure of the Einstein equation in 
the presence of the helical symmetry. In Sec. Ill, we 
briefly describe the gauge conditions which are suited for 
computing gravitational waves from binary neutron stars 
in quasiequilibrium orbits. In Sec. IV, after brief review 
of the conformal flatness approximation and hydrostatic 
equations for a solution of quasiequilibrium states, we 
introduce a linear approximation and derive the equa- 
tions for computation of gravitational waves from the 
quasiequilibrium states. In Sec. V, we numerically com- 
pute gravitational waves from irrotational binary neutron 
stars in quasiequilibrium circular orbits. First, we cali- 
brate our method by comparing the numerical results 
with post Newtonian formulas for gravitational waves 
from two point masses 0J|, adopting weakly gravitat- 
ing binary neutron stars. We will demonstrate that our 
results agree well with post Newtonian analytic formu- 
las Then, gravitational waves from more compact 
binaries are computed to point out the importance of 
tidal deformation and strong general relativistic effects 
on gravitational waves for close binaries. Section VI is 
devoted to a summary and discussion. 

In the following, we use geometrical units in which 
G = c = 1. We adopt spherical polar coordinates; Latin 
indices z, j, k, . . . and Greek indices n,v, . . . take r, 0, ip 
and t, r, 8, <p, respectively. We use the following symbols 
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for a symmetric tensor, A^ — (Aij + Aji)/2 and the 
Kronecker's delta 



II. BASIC EQUATIONS 

We are going to compute gravitational waves from bi- 
nary neutron stars in quasiequilibrium circular orbits us- 
ing an approximate framework of the Einstein equation. 
Before deriving the basic equations for the approxima- 
tion, we describe the full sets of the Einstein equation in 
the presence of a helical Killing vector as 
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(2.1) 



where denotes the orbital angular velocity and l % — 
il(d/d(p)' 1 . The purpose in this section is to clarify the 
structure of the Einstein equation in the helical symmet- 
ric spacetimes. 



A. 3+1 formalism for the Einstein equation 

We adopt the 3+1 formalism for the Einstein equation 
Jl8[ in which the spacetime metric is written as 

ds 2 = g^dx^dx" 

= {-a 2 + f3 J f3 3 )dt 2 + 2/3jdx j dt + lio dx l dx j , (2.2) 

where g M „, a, ft (Pi = Jijft), and 7^ are the 4D met- 
ric, lapse function, shift vector, and 3D spatial metric, 
respectively. Using the unit normal to the 3D spatial 
hypersurface E t , 

n^= (-,-—) and n„ = (-a, 0, 0, 0), (2.3) 
\ a a ) 



jij and the extrinsic curvature K%j are written as 



lij = 9ij + Thill j, 



Kij = -7iS/V fc n ; , 



(2.4) 
(2.5) 



where is the covariant derivative with respect to g^. 
For the following calculation, we define the quantities 

as 



7 = det(7y), 
lij = i>~%j, 



Ay = ip - -JijK 



(2.6) 
(2.7) 

(2.8) 



where tp is a conformal factor and K = Kijj^ . In con- 
trast to the formalism which we use in 3+1 numerical 
simulations Jl9[ , we do not a priori impose the condition 
7 = det(7y ■) = det(?7ij) = rj where rjij is the metric in the 
flat space and rj = r 4 sin 2 9. In the following, the indices 



of variables with a tilde such as Ay, A- 7 , and ft(= ft) 
are raised and lowered in terms of 7y and j lJ . Here, 
Di, and (o)-Di are defined as the covariant derivative with 
respect to 7y, 7y, and 7/y, respectively. 

The Einstein equation is split into the constraint and 
evolution equations. The Hamiltonian and momentum 
constraint equations are 



R - KaK^ +K 2 



D.K'j -D,K = 8-irJj 



or 



A V , = ^R- 2nE^ - (.l, ; .r-' 



-K- 



(2.9) 
(2.10) 



(2.11) 
(2.12) 



where E and Ji are defined from the energy-momentum 
tensor T^ v as 



E = T^n", 



(2.13) 
(2.14) 



Ji = -T^n^Yi ■ 
R and R are the scalar curvatures with respect to 7y and 



7y, and A = D k E) k . The elliptic- type equation (2.11) 
will be used for determining ip. 

The evolution equations for the geometry are 



d t jij= -2aKij + Dify + Djpt, (2.15) 
d t Kij= aRij - DiD ja + aiKKij - 2K ik K j k ) 
+ (D,ft)K H + (Di(3 l )Kij + (DiKij)p 1 
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E-S, 



(2.16) 



where i?y is the Ricci tensor with respect to 7y and 



li k li l T k i- 



(2.17) 



By operating 7^ in Eqs. (2.15) and (2.16), we also have 

^(-aK + D^-tM, (2.18) 
d t K= olK^K 11 -Aa + Ana{E + S k k ) + /^SjJf, (2.19) 

where A = D k D k . To write the evolution equation of 
K in the form of Eq. (2.19), we use the Hamiltonian 
constraint equation (2.9). Using Eqs. ( 2.15Q and (2.18), 
the evolution equation for 7y is described as 



dtlij - ~{dtl)lu 
37 



-2aAy + Di(3j + Djfc - -%jD k ft 



(2.20) 
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B. Einstein equation in helical symmetric spacetime 

In the presence of the helical Killing vector £ M , 7y 
and Kij satisfy -£-^ij = = £$Kij where denotes 
the Lie derivative with respect to £ M . In spherical polar 
coordinates, the relations are explicitly written as 



dtKij= -£ k d k Kij. 



Using Eqs. ( ^2l| ), Eqs. (|I|), ( pH) ), (pJ8| ), and 
are rewritten as 



2aKij= DiLUj + DjUii, 



(2.22) 



0= aiZa - DiDja + ol(KK 13 - 2K lk K k ) 
+ {D 3 J)K H + {DiJ)Kij + JDiK l3 



— Sna 
aK= DiLu\ 



Si 



2 7«( 



— w 
where 



" ' 2 

fe 9kif= aKijK ij -Aa + 47ra(£; + S fe fe ), 



tj fc = /3 fe +£ fe . 



(2.23) 

(2.24) 
(2.25) 



(2.26) 



Equation ( 2.20 ) is also rewritten in the form 



2aAij = DiQj + Dj&i - Aw fc 

= D i i3 j +D j j3 i -pi j Dk0 k 
+i k d k j tJ -^(£ k d k ^, 



(2.27) 



where uji = "fijUj J (Co 1 = lu 1 ), and we have used relation 
dtjij = -i k d k jij- 



Substituting Eq. (|2.27D into Eq. ( p,12[ ), we obtain 
equations for cu l and (3 l as 



Equation (2.29) is solved to determine after we ap- 
propriately specif y the spatial gauge condition for 7^. 
In handling Eq. (2.29), the following relation is useful 
t o eva luate the sum of the fourth and sixth terms in Eq. 
( P9|) : 



Kl l +R\i k =l k (d k f l 



(2.30) 



(2.21) Here, is the Christoffel symbol with respect to 7^. 



The equation for 7^ is derived from Eq. (2.23). For 
the derivation, we first rewrite Rij as 



Rij - Rij + Rij, 



(2.31) 



where Rij is the Ricci tensor with respect to 7^' and 



Rt = 



2 - - 2 
-D i D j rJ J --n fij AiP 

6 - 



+— DiijjDj^ - —jijDhipD^. 



^2 ^ ^2 

Using (o)-Dfc, Rij is written as 
1 



(2.32) 



Rij - 2 



— Aa a t/ijj + ( )Dj(Q}D k h k i + ( jDi( Q -)D k h k j 
- 2( )ACfcj + 2( 0) £) fe (/ fcZ C; ! .y) 



(2.33) 



where Afl at = (o)A (o)D k , and we split 7^ and 7^ as 
rjij + hij and rf 2 + /*■?, respectively. Cy and Ck,ij are 
defined as 



C k j = -^-((o)A/iji + (o)Djhu - ^Dihij^j, 
Gisj = - f (o^Dihji + (tyDjha — (o)A/kA ■ 



(2.34) 



AcDj + -D 3 D k uj k + R jk Cb 



~.k 



and 



D' In 



/"iJj 6 \/~ ~ 2 - ,. 



—aDjK — 16naJj 
o 



A/9j- + ^DjDkj3 k + RjkP k 
+ j jk (Ae k + ^D j D k £ k + R k l l l 



(2.28) 



) 



D l In 



Aft- + Aft " g7« 



+ i k dkli j -^(i k d k 7)jij 
37 



—aDjK — 16iraJj. 



(2.29) 



/ / = ^{ln(V7)} and f,<, 



We note that T\ 
dj{\n(y/ "//?])}. It is also worthy to note that in the linear 
approximation in hij , Li = 7^ L 3 reduces to 



(v)D l h H -^d l (h kl r 1 kl ) +()\th ;i r 



(2.35) 



The second line in Eq. (2.23) is written as 

(DjU> k )K u + (DiOJ k )K kj + u k D k K l3 
= (D^ k )K kl + (D t f3 k )K k] + p k D k Kij 



+ -[Ke k d k ^ + ll3 e k d k K) +td k {^A l3 ). (2.36) 



Substituting Eq. ( 2.27 ) into the last term, we find the 
presence of a term as 

\l k d k {^{t l dihij)). (2.37) 
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Recalling the prese nce o f a term — Afl at /iy/2 in R 



is found that Eq. ( 2.23| ) constitutes a Helmholtz-type 



equation for the nonaxisymmetric wave parts of hij as 



«Aj a 



^ dk) iL^ dl) 

a 



(source) i 



(2.38) 



In the axisymmetric case, the equation for hij changes 
to an elliptic-type equation. This is natural because in 
stationary, axisymmetric spacetime, there do not exist 
gravitational waves. In the nonaxisymmetric case, also, 
the axisymmetric part of hij obeys an elliptic-type equa- 
tion, and hence it is regarded as a nonwave component 

As a consequence of the calculations in this section, it 
appears that ip and (3 l obey elliptic-type equations and 
hence they seem to be nonwave components. However, 
it is not always true. If we would not carefully choose 
gauge conditions, these variables could contain a wave 
component even in the wave zone. To extract gravita- 
tional waves simply from nonaxisymmetric parts of hij, 
it is preferable to suppress wave components in these vari- 
ables with an appropriate choice of gauge conditions. In 
the next section, we propose a gauge condition which 
meets the above demand. 



III. GAUGE CONDITIONS 

In this section, we propose gauge conditions which are 
suited for the computation of gravitational waves emitted 
from quasiequilibrium states. 

As the time slicing, we adopt the maximal slicing con- 
dition as 

K = = d t K. (3.1) 
Then, an elliptic-type equation for a is obtained; 



Aa = 4ira{E + S k k ) + . 



(3.2) 



This equation may be written as 
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aip f, 



A(aV>) = 2na^(E + 2S k k ) + -aifj° A tJ A 11 + -^-R 



(3.3) 



Note that in the case K = 0, it is found from Eq. (2.24) 
that the condition 



D k U) K 







(3.4) 



must be guaranteed in solvi ng E q. (|2.28 ) [or ( 2.29 )]. 
Namely, the solution of Eq. (2.28) in the condition K = 
has to satisfy the relation D k ui k — 0. It is e asy t o 
show that the condition is really guaranteed if Eq. (2.19), 
the Hamiltonian constraint, and the Bianchi identity are 
satisfied. 

We propose spatial gauge conditions for hij in which 



rfihij = 0[(hijf] 



(0) 



D k h 



hi 



(0) 



D*Tn 



and 

'V> 6 ' 



h ki = 0[(hijf 



(3.5) 



where on the right-hand side of these equations, we allow 
adding certain nonlinear terms of h^. For simplicity, 
we consider here the case in which they are vanishing. 
Namely, we adopt a transverse and tracefree condition 
for ip 6 hij/a. In this case, 



7 = r / {l + 0[(^) 2 ]}. 



(3.6) 



There are two merits in choosi ng thi s gauge condition. 
The first one is that using Eq. (2.35), we can derive a 
relation in this gauge as 



Lj + D l In 



v> 6 



-e k d k 



(o)Dj\n 



V> 6 



0[{hi 



(3.7) 



Thus, the equation for determining f3 l [Eq. (2.29)] does 



not contain linear terms of hij except for coupling terms 
between /3 l and hij and between d v \di ln(ip 6 /a)] and hij. 
Since the magnitude of these coupling terms and non- 
linear terms of hij is much smaller than that of leading 
order terms such as Afl at /3i and I6ira Ji, we can consider 
that effects due to hij are insignificant in the solution 
of /3*. If information on gravitational waves is mainly 
carried by hij, not by other metric components, the so- 
lution of the equation for (3 l is not contaminated much 
by the wave components and it is mainly composed of 
a nonwave component in the wave zone. As a result of 
this fact, it is allowed to regard (3 k in the wave zone as a 
nonwave component. 

In the maximal slicing condition K — 0, the following 
relation holds: 



= -^d k [^{l k 



+ (3 k )}. (3.8) 



Since the right-hand side of this equation is weakly de- 
pendent on hij and mainly composed of nonwave com- 
ponents, we may also regard ip in the wave zone as a 
nonwave component. 

The second merit appears in the equation for h^ , which 
is written as 



lat--(^dfc) — (^) h i3 
a a J 

+ ^)Dii[hj )m D k 
2 {-(o) D i C k] + (a)D k {f kl Ciij) - CljC^ 
+ C\jC k kl + R% } - ^DiDja - ^A lk A k J 
+ ^\2^A k{l bj{p k + p k D k (^A l3 )} 
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1 

a 



Q 



+ -l k d k { ( Dipj + D 3 (3 t - =-%D n f3 n 



37 



8n[2S ij+lij (E-S k k )}, 



(3.9) 



where we use the condition K = 0. Osn the left-hand 
side, only linear terms in hij are collected, and on the 
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right-hand side, the nonlinear terms are located. (Note 
' - 0{hij) and C|, = 0[{hij) 2 }.) In the lin- 



that C& 



ear order in hij, Eq. (3J)) is regarded as a Helmholtz- 



type equation in a curved spacetime for nonaxisymmetric 
parts of h^. As a result, we can clarify that the nonax- 
isymmetric parts of hij are wave components in the wave 
zone. This fact is helpful in specifying the boundary con- 
dition in the wave zone. 

Since both wave and nonwave components are in- 
cluded, it is not trivial how to impose outer boundary 
conditions for hij . A solution to this problem is to use a 
spectrum decomposition method in which we expand h%j 
as 



h[^ exp(im< y £>), 



(3.10) 



and solve each m mode separately. As already clarified, 
h\j is a nonwave component and h\j (m 7^ 0) is a wave 
component. Thus, we can impose the outer boundary 
condition for both components correctly. 

Before closing this section, the following fact should be 
pointed out. For computation of quasiequilibrium states 
in the presence of the helical Killing vector, the minimal 
distortion gauge |2l| in which 



AW> 4 7 1/3 9t¥ J ') = 



(3.11) 



is not available. In this gauge, we fix the gauge condition 
for dtjij, but do not specify any gauge condition for 7^; 
i.e., an initial gauge condition at t = is not specified. 
To obtain a quasiequilibrium state, on the other hand, we 
have to fix the gauge condition initially, and as a result, 
throughout the whole evolution, the gauge condition is 
fixed because of the presence of the helical Killing vec- 
tor. This is the reason that we cannot use the minimal 
distortion gauge in the helical symmetric spacetimes. 



IV. FORMULATION FOR COMPUTATION OF 
GRAVITATIONAL WAVES 

A. Equations for background quasiequilibrium 
neutron stars 

Instead of solving the full equations derived above, in 
this paper, we adopt an approximate method for the 
computation of gravitational waves from binary neutron 
stars in quasiequilibrium states. First, we compute the 
quasiequilibrium states of binary neutron stars in the 



framework of the so-called conformal flatness approxima- 
tion neglecting hij [p2 16,17[|. Then the basic equations 
for the gravitational field are 



A flat (aV>) = 2naiP 5 (E + 2S k k ) + A t] A^ , 



V> 5 7 



Aij A 1 - 1 , 



A fl at/3, + ^ 0) D m D k p k + ,, I'^-'i^ ) 
= l6iraJj, 

where 
A 



(o)Di/3j + (o)Dj/3i - -r)ij(o)Dkl3 , 



1 



(4.1) 
(4.2) 

(4.3) 

(4.4) 
(4.5) 



and we set K = 0. The spatial gauge condition (3.5) is 



automatically satisfied since we assume hi 



0. 



In the far zone, these gravitational fields behave as 
M 



a - 



1 h y J ai m {r)Yi m , 

l>2,m 
l>2,m 

hm{r)(Y lm , 0,0) 

l>l.m 



1 



(4.6) 
(4.7) 



T L7TL 

+c lm (r) (0, 8 v Y lrn / sin 6, -d e Y lrn sin 9)] , (4.8) 

where M denotes the Arnowitt-Deser-Misner (ADM) 
mass of the system, and Yi m (9,(f) is the spherical har- 
monic function. We implicitly assume that the real part 
of Yi m is taken. The asymptotic behaviors of ai m , ipim, 
ai m , hm> and c Zm at r — > 00 are 



aim 

Iplm 
C-lm 



-l-l 



-l-l 



-1-2 



(4.9) 



The coefficient of the monopole part of a should be — M 
for quasiequilibrium states in the conformal flatness ap- 
proximation pq ] . This relation is equivalent to the scalar 
virial relation so that it can be used for checking numer- 
ical accuracy [see Eq. (5.7)]. 

We adopt the energy-momentum tensor for the perfect 
fluid in the form 



Tp U = {p + pe + P)u,j,u v + Pg^ v , 



(4.10) 



where p, e, P, and denote the rest mass density, spe- 
cific internal energy, pressure, and four-velocity, respec- 
tively. We adopt polytropic equations of state as 



G 



P : 



up 



(4.11) 



where n is a polytropic constant, T = 1 + 1/n, and n a 
polytropic i ndex . Using the first law of thermodynamics 
with Eq. ( 4.11 ), e is written as nP/p. The assump- 
tion that k is constant during the late inspiraling phase 
is reasonable because the timescale of orbital evolution 
for binary neutron stars due to the radiation reaction of 
gravitational waves is much shorter than the heating and 
cooling timescales of neutron stars. In this paper, we 
adopt n — 1 as a reasonable qualitative approximation 
to a moderately stiff, nuclear equation of state. 

Since the timescale of viscous angular momentum 
transfer in the neutron star is much longer than the evo- 
lution timescale associated with gravitational radiation, 
the vorticity of the system conserves in the late inspi- 
raling phase of binary neutron stars p3]. Furthermore, 
the orbital period just before the merger is about 2 ms 
which is much shorter than the spin period of most of 
neutron stars. These imply that even if the spin of neu- 
tron stars would exist at a distant orbit and would con- 
serve throughout the subsequent evolution, it is negligi- 
ble at close orbits for most of neutron stars of the spin 
rotational period longer than ~ 10 ms. Thus, it is rea- 
sonable to assume that the velocity field of neutron stars 
in binary just before the merger is irrotational. 

In the irrotational fluid, the spatial component of 
is written as 



h 



(4.12) 



where h = 1 + e + P/ p and $ denotes the velocity po- 
tential. Then, the continuity equation is rewritten to an 
elliptic-type equation for $ as 



Diipah^D*®) - D i [pau t {e i + = 0. 



(4.13) 



In the presence of the helical Killing vector, the relativis- 
tic Euler equation for irrotational fluids can be integrated 
to give a first integral of the Euler equation as J25| 



+ hukV k = const, 

u 



(4.14) 



where V k = u k /u l - £ k . Thus, Eqs. fl4.13| ) and fl4.14| ) 
constitute the basic equations for hydrostatics. 



B. Equation for hi 



After we obtain the quasi equ ili briu m state s solvi ng th e 
coupled equations of Eqs. (|jHO), fl4.13| ), and (M4 



the wave equation for hij [Eq. ( |3.9| )] is solved up to linear 
order in hij in the background spacetime of the quasiequi- 
librium states. Without linearization, nonlinear terms of 
hij cause a problem in integrating the equation for hij in 
the wave zone because standing gravitational waves exist 
in the wave zone in the helical symmetric spacetimes and 



as a result the nonlinear terms of hij fall off slowly as r~ 2 . 
In a real spacetime, the helical symmetry is violated be- 
cause of the existence of a radiation reaction to the orbits. 
This implies that the existence of the standing wave and 
the associated problem are unphysical. Thus, we could 
mention that linearization is a prescription to exclude an 
unphysical pathology associated with the existence of the 
standing wave. 

In the absence of nonlinear terms of gravitational 
waves, we cannot take into account the nonlinear memory 
effect [^6j. However, as shown in p6fl , this effect builds 
up over a long-term inspiraling timescale, and as a result, 
it only slightly modifies the wave amplitude and luminos- 
ity of gravitational waves at a given moment. Thus, it is 
unlikely that its neglect significantly affects the following 
results. 

In addition to a linear approximation with respect to 
hij, we carry out a further approximation, neglecting 
terms of tiny contributions such as coupling terms be- 
tween (3 k and hij and between T M „ and h^. We have 
found that the magnitude of these terms is much smaller 
than the leading order terms and its contribution to the 
amplitude of gravitational waves appears to be much 
smaller than the typical numerical error in this paper 
of ~ 1%. We only include coupling terms between hij 
and spherical parts of a and ip since they yield the tail 
effect for gravitational waves which significantly modi- 
fies the amplitudes of gravitational waves |27|| . We also 
neglect the perturbed terms of ip, a, and /3' associated 
with h^ since they do not contain information of gravita- 
tional waves in the wave zone under the gauge conditions 
adopted in this paper |p8| . With these simplifications, 
the numerical procedure for a solution of hij is greatly 
simplified. 

As a consequence of the above approximation, we ob- 
tain the wave equation of hij as 



A 



flat 



^{i k d k f 



2Rt i --DiD J a-4i> 4 A lk A k 1 



a 



^{2^Ak {m D. }) p k + p k {0) D k {^A l3 )) 

V> 4 



i k d k 



-{L(3)i 



8Tr[2Sij+<tp\j(E-S k k )} 



QE 



8R y JDiDja 
13 \ a 



(4.15) 



where [• • -]qe is calculated by substituting the geomet- 
ric and matter variables of quasiequilibrium states: In 
the following, we denote it as S^ E . Here SRfj and 
6(DiDja/a) denote coupling terms between linear terms 
of h^ and ipo or ao in Rf- and DiDjCt/a, and 4>o and oiq 
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denote the spherical part of ip and a which are computed 
by performing the surface integral over a sphere of fixed 
radial coordinates as 



Qo(r) = ^-<f QdS, 

J r— const. 



(4.16) 



where dS = sin 6d9d<p. Note that in the present formu- 
lation, the spatial gauge condition is 



rf 3 hij = and 



{0) D K h H + \ (o)-D In 



(5) 



h ki = 0. 



(4.17) 



We neglect coupling terms of hij with a and ip except for 
with «o and ipo since their order of magnitude is as small 
as that of coupling terms between hij and /3 1 



In Eqs. Q4.15| ) and fl4.17| ), V'o/^o and ^/"o are dif- 
ferent from the spherical part of 4> 6 /a and ip A /a 2 in the 
near zone, although in the wave zone they a re al most 
identical. This implies that in deriving Eqs. (4.15) and 
(4.17), certain ambiguity remains. However, in numer- 
ical computations, we have found that the difference of 
the numerical results between two formulations is much 
smaller than the typical numerical error. For this reason, 
we fix the formulation and gauge condition in the form 
of Eqs. ( |4~15| ) and ( p^ ). 

As mentioned in Sec. I, the procedure for the compu- 
tation of gravitational waves adopted here is quite simi- 
lar to that for obtaining gravitational waves in the post 
Newtonian approximation H|J: In the post Newtonian 
work, one first determines an equilibrium circular orbit 
from post Newtonian equations of motion, neglecting the 
dissipation terms due to gravitational radiation. Then, 
one substitutes the spacetime metric and matter fields 
into the source term for a wave equation of gravitational 
waves. In this case, no term with regard to gravitational 
waves and the radiation reaction metric is involved in the 
source term of the wave equation. We may explain that 
we here follow this procedure. 

The main difference between our present method and 
post Newtonian calculations is that we fully include a 
general relativistic effect in ip, a, and (3 k for quasiequi- 
librium states without any approximation, and that we 
take into account the effect of tidal deformation of each 
neutron star. 

In the post Newtonian approximation, hq is p resent 
from second post Newtonian (2PN) order |pO| , pl| . This 
implies that quasiequilibrium states obtained in the con- 
formal flatness approximation and gravitational waves 
computed in their background spacetimes contain an er- 
ror of 2PN order from the viewpoint of the post Newto- 
nian approximation. To obtain quasiequilibrium states 
and associated gravitational waveforms for better post 
Newtonian accuracy, it is necessary to take into account 
h^. Here, we emphasize that our method does not re- 
strict the zeroth-order solution of the three-metric in con- 
formally flat form. Even if quasiequilibrium states are 



constructed in a formalism with hij, we can compute 
gravitational waves in the same framework. In this pa- 
per, we adopt the conformal flatness approximation sim- 
ply because of a pragmatic reason as mentioned in Sec. 
I. With a modified formalism and a new numerical code 
taking into account hij , it would be possible to improve 
the accuracy of quasiequilibrium states appropriately in 
the present framework (see discussion in Sec. VI). 



C. Basic equations for computation of hij 



Since hij in Eq. ( 4.15| ) couples only with functions of 
r, we decompose it using the spherical harmonic function 
Yi m (6, if) as 



l,m 
+ Blm 

+ ]>> 2 F fa 

1, 171 

l,m 



Yl m 
,2i 





* -r 2 Y lm /2 

* * -r 2 sm 2 0Yi m /2 

III 

* 

* * 





* Wl„ 

* * 





sin 2 9Wi„ 

d v Y lm /sin9 -d g Y hn smt 

* 

* * 



* -X hn /sm9 W hn sm9 

* * sin 9Xi 

ill 



(4.18) 



where * denotes the relations of symmetry. Note th at th e 
trace- free condition for hij is used in defining Eq. ( 4.18 ). 
Here, A lm , B lm , C; m , D tm , and F lm are functions of r, 
and 



W lm = 
Xlm = 2d v 



(def - cot9d e - -^-(c^) 2 Y lm , (4.19) 
sin 9 



de — cot 9 



Y,, 



(4.20) 



Using Eq. ( |4.18] ), the equations of the spatial gauge 
condition are explicitly written as 



dAi m 3 

—j r - Air, 

ar r 



M — 5 r Ai m 



dr \a J 



dB, 



- r - + -B lm -^-\,F, 
ar r 2 



l r lm 



dC, 



lm 



dr 



r 



dr \a J 

dr \a J 



= 0, 

6' 



= 0, (4.21) 

(4.22) 
= 0, (4.23) 



where A/ = 1(1 + 1) and A; = A/ — 2. From these equations, 
we find that the following relations have to be satisfied 
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in this gauge condition: (1) for I = 0, A\ m cx ao/r 3 -0Q 
with B lm = F lm = C'i m = D im = 0; (2) for I = 1, 
Ai m cx a /r 2 ip$ (B lm cx a /rip$) or A hn cx a /r 4 ip$ 
(B lm oc a /r 3 ip§) with F hn = 0; (3) for I = 1, C im cx 
oio/r ipQ with -D; m = 0. The behavior of Ai m ,Bi m , and 
C; m for ^ = and 1 is regular for r -> oo, but not for 
r — * 0. This implies that they should vanish for I = 
and 1, and modes only of I > 2 should be nonzero. Thus, 
nonwave components in hij of I < 1 can be erased in the 
present gauge condition. 

For I > 2, Bj m , jFj m ) and i?; m can be calculated from 
Ai m and Ci m because of our choice of the spatial gauge 
condition. This implies that we only need to solve the 
equations for Ai m and C; TO , which are derived as (cf. 
Appendix A) 



dr 2 



2 ± In 

dr 



(5) 



d_ 

dr 



A; — 6 



^0 2r>2 

—3- mi 2 



2 {dY 2 ln (a 



A, 



d In a?o dA 



I in 



dr 



dr 

dAi„ 
dr 



t dla 

dr 
6Ai m 



2A, 



' dr 



In 



(^8 



(4.24) 



dr 2 { r dr \ ao J J dr r 2 



2 d W / 
r dr 



41 — \nipo 
ar 



+4 



d 



1 

qC 



dr 

(is [ s: 



\ni/j Q 



dr 



In «o 



Clr, 



smf 



(4.25) 



where Y,^ denot es th e complex conjugate of Yi mi and 
Eqs. ( 4.21 ) and ( 4.23 ) are used to erase Bi mi F\ m and 
Di m in these equations. For the case m = 0, these equa- 
tions are elliptic- type equations for Ai m and C; m , imply- 
ing that they are not gravitational waves. 

In this paper, we consider the binaries of two identical 
neutron stars. Then, the system has w- rotation symme- 
try. In this case, Ai m of even I and even m and Ci m of 
odd / and even m are nonzero, and other components are 
zero. 

S{m and Sim of to 7^ behave as 0(/~' - 1 ) fo r r — > 00 
because of the presence of a term [cf. Eq. (O)] 



a \ a 



(4.26) 



For I — 2 and 3, the falloff of this term is so slow that it 
could becom e a so urce of numerical errors in integrating 
Eqs. (4.24) and (4.25) for the computation of gravita- 
tional waves in the wave zone. Furthermore, Eq. ( 4.26 ) 
gives a main contribution for solutions of Ai m and C/ m 
in the wave zone; namely, we need to carefully estimate 
the contribution from this term for an accurate compu- 
tation of gravitational waves. To resolve this problem, 
we transform the variables from Ai m and C; m to new 
variables as 



.4 



lm 



A 



1 



lm 



imQ ./ r = cons t. 



dS(LP) rr Yi m , 



(4.27) 



Cim — Cim + 



1 



dS 



l J r— const. 



r) Y* 



, (4.28) 



and rewrite Eqs. (4.24) and (4.25) in terms of Ai m and 
Cim- With this procedure, the source terms of the wave 
equations for Ai m and C/ m fall off as 0(r~ l ~ 3 ), so that it 
becomes feasible to accurately integrate the wave equa- 
tions without technical difficulty. 



D. Boundary conditions 

Ordinary differential wave equations for Ai m and Ci rn 
with 2 < / < 6 and 2 < |m| < 6 are solved, imposing 
boundary conditions at r = as 



dAim dCi 

T: 



dr 



dr 







(4.29) 



and at a sufficiently large radius r = r max 3> A = 
27r(mf2) -1 as 



d(r 3 Aim) 
dr* 

djrCim) _ 
dr* 



= imflr 3 Aim, 
= imflrCim, 



where r* denotes a tortoise coordinate defined by 

r* = fdr^l. 
J a 

Here, we assume the asymptotic behaviors 
exp(imf2r*) 



Aim — * Ca- 



Cim —> Cc- 



exp(im51r*) 



(4.30) 
(4.31) 

(4.32) 

(4.33) 
(4.34) 



where Ca and Cc are constants. 

Note that for obtaining an "equilibrium" state in which 
no energy is lost from the system, we should adopt the 
ingoing-outgoing wave boundary condition for keeping an 
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orbit. However, the purpose here is to compute realis- 
t ic, o utgoin g gra vitational waves, so that we adopt Eqs. 
(4.33) and (4.34) as the outer boundar y con ditions. 

For m — 0, the falloff of the term (|4.26|) is not very 
slow, so that we do not have to change variables. Elliptic- 
type ordinary differential equations (ODEs) for Ai m and 
Cim are solved, imposing the boundary conditions at r — 
as 



dr 

d(C hn /r l ) 



= 0, 



dr 



= 



and the outer boundary conditions as 



Aw, 



-1-2 



-l-l 



(4.35) 
(4.36) 



(4.37) 
(4.38) 



Note that these outer boundary conditions are deter- 
mined from the asymptotic behavior of their source terms 
[cf. Sf 3 E and Eq. Q]. 

Since the wave equations are ODEs, it is easy to take 
a sufficiently large number of grid points up to a dis- 
tant wave zone in current computational resources. If 
the outer boundary conditions are imposed in the dis- 
tant wave zone, the above simple boundary conditions, 
without including higher order terms in 1/r, are accept- 
able. Also, ODEs can be solved with a very high accuracy 
in current computational resources. Thus, the numerical 
accuracy for gravitational waveforms computed below is 
limited by the accuracy of quasiequilibrium states ob- 
tained in the first step (i.e., the source terms of the wave 
equations limit the accuracy). 



E. Formulas for gravitational wave amplitude and 
luminosity 

In the distant wave zone, + and x modes of gravita- 
tional waves, h + and h x , are defined as [M 



1 

2r~ 2 



sin 2 6 



h x = 



and, thus, 



2<Z<6 



Y, [FlmWl m -D h 



sin I 



h x = ^ ( Flm- 



2</<6 



Xl, 



DlmWlr, 



(4.39) 
(4.40) 

(4.41) 
(4.42) 



In the distant wave zone, Fi m and Di m can be obtained 
from Ai m and Ci m as 



F,,. 



{mn) 2 Ai m r 2 

Xi~Xi 
imQCim 



X 



(4.43) 
(4.44) 



For the latter, we write h + in the wave zone as 



h+ = - 



#22(1 + u 2 ) cos(24<) + H 32 (2u 2 - 1) 008(2*) 



+H i2 (7u 4 - 6u 2 + 1) cos(2*) 
+#44(1 - u 4 ) cos(4*) 
+ J ff 5 2(12u 4 - Uu 2 + 1) cos(2*) 
+iT 54 (4u 2 - 1)(1 - u 2 ) cos(4#) 
+ J ff 62 (495u 6 - 735u 4 + 289u 2 - 17) cos(2*) 
+H 6A [33u 4 - 10u 2 + 1)(1 - u 2 ) cos(4#) 

+H 66 (u 2 + 1)(1 -u 2 ) 2 cos(6*) , (4.45) 

where * = <p — fit, D is the distance from a source to 
an observer, and Hi m denotes the amplitude for each 
multipole component (/,m). Here, we assume that the 
mass centers for two stars are located along x-axis at 
t = 0. The gravitational wave luminosity is computed 
from H 

§ = ^-1 dS(h 2 + + hi) 

dt lD7T ,/ r .„ 



D 2 Xi~Xi 
16tt 



(mrt) 2 (\Fi m \ 2 + \Di m \ 2 ), (4.46) 



where h-> 



2<l<6 

dh +iX /dt. 



V. NUMERICAL COMPUTATION 

A. Numerical method and definition of quantities 

1. Computation of zeroth- order solutions : Quasiequilibrium 
sequence of binary neutron stars 

Following previous works @0, we define the coor- 
dinate length of semimajor axis Ro and half of orbital 
separation d for a binary of identical neutron stars as 



Ro 



d = 



Rout R'm 



Rout + R in 



(5.1) 
(5.2) 



where R m and i? ut denote coordinate distances from the 
mass center of the system (origin) to the inner and outer 
edges of the stars along the major axis. To specify a 
model along a quasiequilibrium sequence, we in addition 
define a nondimensional separation as 
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d = — 



(5.3) 



At d = 1, the surfaces of two stars contact and at d — > oo, 
the separation of two stars is infinite. In the case n = 1, 
the sequences of binaries terminate at d = o? m in — 1-25 
for which the cusps (i.e., Lagrangian points) appear at 
the inner edges of neutron stars jl7|| . Also it is found 
that for d > 2, the tidal effect is not very important. 
Thus, we perform a computation for 1.25 < d < 3. 

In using the polytropic equations of state (with the 
geometrical units c = G = 1), all quantities can be nor- 
malized using k as nondimensional as 



M = Mk" /2 , J = Jn n , 



(5.4) 



where M, J, and i? denote the total ADM mass, total 
angular momentum, and a circumferential radius. Hence, 
in the following, we use the unit with k = 1. For later 
convenience, wc also define several masses as follows: 

Mq : the rest mass of a spherical star in isolation, 
M g : the ADM mass of a spherical star in isolation, 
M t = 2M g , 

M : the total ADM mass of a binary system. 

Here M is obtained by computing the volume integral of 
the right-hand side of Eq. (4.2). Note that M is not equal 
to M t in the presence of the binding energy between two 
stars. 

The binding energy of one star in isolation and the 
total binding energy of the system is defined as 



E b 
E t 



M g - M , 
M - 2M . 



(5.5) 
(5.6) 



The energy and angular momentum are monotonically 
decreasing functions of d(> d m i n ) for n = 1 irrespec- 
tive of the compactness of each star. 

Quasiequilibrium states in the framework of the con- 
formal flatness approximation are computed using the 
method developed by Uryu and Eriguchi Q . We adopt 
a spherical polar coordinate (r, 9, ip) in s olvi n g ba sic equa- 
tions for gravitational fields [cf. Eqs. ( 4.1 )— ( 4.3 )] . Here, 
the coordinate origin is located at the mass center of the 
binary. Since we consider binaries of identical stars, the 
equations are numerically solved for an octant region as 
< r < 100.Ro and < 0,ip < tt/2. We typically take 
uniform grids of 51 grid points for 9 and ip. For the ra- 
dial direction, we adopt a nonuniform grid and the typical 
grid setting is as follows: For < r < 5i?o, we take 201 
grid points uniformly (i.e., grid spacing Ar = 0.025i?o)- 
On the other hand, for 5i? < r < 100i? , we take 240 
nonuniform grids, i.e., in total 441 grid points for r. A 
fourth-order accurate method is used for finite differenc- 
ing of 9 and p directions and a second-order accurate one 



is used for r direction. Hydrostatic equations are solved 
using the so-called body- fitted coordinates (r',9', ip') (l(J 
which cover the neutron star interior as < r' < Ro, 
< 9' < tt/2, and < ip' < n, respectively. We adopt 
a uniform grid spacing for these coordinates with typi- 
cal grid sizes of 41 for r', 33 for 9', and 21-25 for ip' . 
A second-order accurate finite differencing is applied for 
solving the hydrostatic equations. 

Using this numerical scheme, we compute several se- 
quences, fixing the rest mass M and changing the binary 
separation d. Such sequences are considered to be evo- 
lution sequences of binary neutron stars as a result of 
gravitational wave emission. We characterize each se- 
quence by the compactness which is defined as the ratio 
of the gravitational mass M g to the circumferential ra- 
dius R c of one star at infinite separation. Hereafter, we 
denote it as (M/R)^ [cf. Table I for relations between 
M g , Mo, and (M/R)^]. Computations are performed 
for small compactness {M/R)^ — 0.05 for calibration as 
well as for realistic compactness as (M/R)^ = 0.14 and 
0.19. Relevant quantities of each sequence are tabulated 
in Tables II and III. 

Convergence of a numerical solution with increasing 
grid numbers has been checked to be well achieved. Some 
of the results are shown in jl6| so that we do not touch on 
this subject in this paper. In addition to the convergence 
test, we also check whether a virial relation is satisfied in 
numerical solutions: In the framework of the conformal 
flatness approximation, the virial relation can be written 
in the form fl23| 



VE 



2a^S k k 



d 6 x = 0. 



(5.7) 



As mentioned above, this relation is equivalent to that 
where the monopole part of a is equal to — M. Since 
this identity is not trivially satisfied in numerical solu- 
tions, violation of this relation can be used to estimate 
the magnitude of numerical error. The nondimensional 
quantity VE/M is tabulated in Tables II and III, which 
are typically of O(10 -5 ). We consider that this is sat- 
isfactorily small so that the quasiequilibrium states can 
be used as zeroth-order solutions for the computation of 
gravitational waves. 

The computations in this paper can be carried out even 
without supercomputers. We use modern workstations 
in which the typical memory and computational speed 
are 1 Gbyte and several 100 Mflops. Numerical solutions 
of quasiequilibrium states are obtained after 350 - 650 
iteration processes. For one iteration, it takes about 50 
sec for a single Dec Alpha 667MHz processor so that 
about 7-10 hours are taken for computation of one model. 
With these computational resources, the computation in 
this paper has been done in one month. 
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2. Computation of ODEs for hij 

For solving one-dimensional wave equations for Aj m 
and Cim a uniform grid with the grid spacing Ar and 10 5 
grid points is used. The outer boundary is located in the 
distant wave zone as ~ 80d~ 3/2 A in this setting. Thi s 
make s the simple outgoing boundary conditions ( 4.30| ) 
and ( 4.31 ) appropriate (see discussion below). To ob- 
tain Sfajr) and S^ m {r) in every grid point, appropriate 
interpolation and extrapolation are used. The extrapo- 
lation for r > 100i?o is performed taking into account 
t he asymptotic behavior for a, (3 k , and ip shown in Eq. 
(4.9). The equations for Ai m and Ci m are solved by a 
second-order finite-differencing scheme jointly used with 
a matrix inversion for a tridiagonal matrix |32|] . One- 
dimensional elliptic-type equations for Aio and Cw are 
solved in the same grid setting, only changing the outer 
boundary conditions. These numerical computations can 
be performed in a few minutes using the same worksta- 
tion described above. 



B. Calibration of gravitational wave amplitude and 
luminosity 

1. Convergence test 

Convergence tests for the gravitational wave amplitude 
have been performed, changing the resolution for the 
computation of quasiequilibrium states for every com- 
pactness. As mentioned above, the error associated with 
the method for integrating the one-dimensional wave 
equation is negligible. Since the source terms of the wave 
equation are composed of quasiequilibrium solutions, the 
resolution for the quasiequilibrium affects the numerical 
results on gravitational waves. To find the magnitude 
of the numerical error, the grid size is varied from 51 
to 41 and 61 for 9 and ip and from 441 to 221 and 331 
for r. It is found that varying the angular grid resolu- 
tion very weakly affects the numerical results within this 
range; the convergence of the wave amplitude is achieved 
within ~ 0.1% error. The effect of the varying radial 
grid size is relatively large, but we find that with a typ- 
ical grid size of 441, the numerical error for the wave 
amplitude is < 1% for (I, m) = (2, 2) and < 2% for (3,2), 
(4,2), and (4,4). Since the amplitude of the (2,2) mode is 
underestimated by < 1%, in the following, the total am- 
plitude and luminosity of gravitational waves are likely 
to be underestimated by < 1% and < 2%, respectively. 



2. Comparison between numerical results and post 
Newtonian formulas for a weakly gravitating binary 

Before a detailed analysis on gravitational waves from 
compact binary neutron stars, we carry out a calibration 
of our method and our numerical code by comparing the 



numerical results with the post Newtonian formulas for 
a binary of small compactness (M/R)^. For calibration 
here, we adopt (M/R)^ = 0.05 (cf. Table I for M g and 
Mo, and Table II for the quasiequilibrium sequence). 

We compare the numerical results with post Newto- 
nian formulas of gravitational waveforms for a binary of 
two point masses in circular orbits. Defining an orbital 
velocity as v = (MjO) 1 / 3 , the post Newtonian waveform 
from the two point masses orbiting in the equatorial plane 
is decomposed in the form P,p3[ 



h + = 



2i]Mv 



D 



J? 22 (l + u 2 )cos(2*) 

H 32 (2u 2 - l)cos(2#) 

#42 (7u 4 - 6u 2 + 1) cos(2*) 

H M (1 -u 4 )cos(4*) 
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H u {Au 2 - 1)(1 - u 2 ) cos(4*) 

# 62 (495u 6 - 735?/ + 289/ - 17) cos(2*) 

i?64(33u 4 - 10u 2 + 1)(1 - u 2 ) cos(4*) 



+ j H' 66 (m 2 + 1)(1-m 2 ) 2 cos(6*) 



(5.8) 



where u = cos 9, rj denotes the ratio of the reduced mass 
to M t which is 1/4 for equal-mass binaries, and 
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(5.9) 



Here the post Newtonian order of the modes with / > 7 
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is higher than the third post Newtonian order, and such 
modes have not been published @|. We note that in H22, 
H32, H42, and H44, we include the effect of their tail 
terms of second and half post Newtonian (2.5PN) or- 
der which could give a non-negligible contribution to the 
wave amplitudes. These terms have not been explicitly 
presented in any papers such as J2j]^], but those for H32, 
H42, and H44 may be guessed from black hole pertur- 
bation theory ||] and that for H 2 2 are computed with 
help o f the 2. 5PN gravitational wave luminosity [see Eq. 
( |5.10D ] and Eq. Q4.46Q . The x mode can be written in 
the same way in terms of Hi m , simply changing the de- 
pendence of the angular functions. Hence, we hereafter 
pay attention only to Hi m in comparison. 

We also compare the numerical gravitational wave lu- 
minosity with the 2.5PN formula p,p4[ 
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~dt 
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For r) = 1/4, the first post Newtonian (1PN), 2PN, and 
2.5PN coefficients are -373/84, -59/567 and -373tt/21, 
respectively. Since the 2PN coefficient is by chance much 
smaller than others, the 2PN formula is not different from 
1.5PN formula very much for equal-mass binaries. 

Before we perform the comparison between numeri- 
cal results and post Newtonian gravitational waves, we 
summarize possible sources of the discrepancy between 
two results. One is associated with the conformal flat- 
ness approximation adopted in obtaining quasiequilib- 
rium states. In this approximation, we discard some 
terms which are as large as a 2PN term from viewpoint 
of the post Newtonian approximation. As a result, the 
magnitude of the difference between two results could be 
of 0(v 4 ). The second source is purely a numerical error 
associated with the finite differencing. The magnitude of 
this error will be assessed in the next subsection. The 
third one is associated with the post Newtonian formu- 
las in which higher order corrections are neglected. This 
could be significant for binaries of large compactness. In 
the following, we will often refer to these sources of dis- 
crepancy. 

Calibration for the gravitational wave amplitude 

In Fig. 1, we show the relative difference of Hi rn to 
2rjMv 2 Hi m as a function of v 2 . Here, the relative differ- 
ence is defined as 



RE 



I III 



2r 1 Mv 2 Hi r , 



(5.11) 



The data points are taken at d ~ 1.3, 1.4, 1.6, 1.8, 2.0, 
2.2, 2.6, and 3.0, and v 2 is roughly equal to (R/M) OQ /d. 



We do not consider (5, 2) and (6, 2) modes because their 
magnitude is much smaller than that of the (2, 2) mode. 

We plot three curves for the (2, 2) mode; one curve 
(dotted line) is plotted using the 2PN formula of H22 
shown in Eq. (|5~9|), the second one (solid line) using 
the 1.5PN formula neglecting the 2PN and 2.5PN terms, 
and the third one (thin solid line) is using the Newtonian 
formula [labeled by (2.2)iV]. By comparing the relative 
errors for (2, 2) modes with three post Newtonian formu- 
las, it is found that the post Newtonian corrections up to 
1.5PN order give a certain contribution by ~ 3% of the 
leading order Newtonian term even at d ~ 3 (v 2 ~ 0.017) 
but that 2PN effects are not very important for small 
compactness (M/R)^ = 0.05. It is reasonable to expect 
that post Newtonian correction terms higher than 2PN 
order beyond the leading terms are also unimportant for 
other modes wit h this compactness. This indicates that 
Hi m in Eq. ( |5.9D contains sufficient correction terms for 
I = 2, 3, and 4. On the other hand, the absence of post 
Newtonian correction terms beyond the leading term in 
Hi m f° r I — 5 and 6 would cause an error of a certain 
magnitude (see below). 

The result presented here also indicates that system- 
atic error associated with the conformal flatness approx- 
imation for background binary solutions, in which we 
neglect h t j of 2PN order, is likely to be irrelevant for 
(M/i?)oo = 0.05. 

For I — 2,3, and 4 modes at sufficiently large separa- 
tion as d ~ 3 (v 2 ~ 0.017) in which post Newtonian cor- 
rections and tidal deformation effects become unimpor- 
tant, the relative errors converge to constants as shown 
in Fig. 1. These constants can be regarded as a nu- 
merical error because they should be zero for sufficiently 
distant orbits. Thus, we can estimate that the magni- 
tude of the numerical error is < 1% for (I, in) = (2,2), 
~ 3% for (3,2), and - 1 - 2% for I = 4. These results 
are consistent with those for convergence tests. 

For I = 5 and 6, the post Newtonian formulas we use 
in this paper are not good enough as a theoretical pre- 
diction. Observing the results for I = m = 2 in Fig. 
1, th e post Newtonian formulas for I — 5 and 6 in Eq. 
fl5.9| ) overestimate the true value of the wave amplitude 
by - 3% at d ~ 3 {v 2 ~ 0.017) because of the lack of 
correction terms of 0(v 2 ) and 0(v 3 ) to the leading term. 
Taking into account this correction, we may expect that 
the numerical errors are ~ 4% for I = 5 and ~ 2% for 
1 = 6. These results indicate that our method can yield 
fairly accurate waveforms of gravitational waves even for 
higher multipole modes. 

With decreasing the orbital separation, the ratio of the 
numerical to post Newtonian amplitude becomes higher 
and higher irrespective of (l,m). This amplification is 
due to the tidal deformation of each star |pq] . For the 
(2, 2) mode, the amplification factor is not very large, 
i.e., ~ 2%, even at d = 1.3 (v 2 ~ 0.035). However, 
for higher multipole modes, the amplification factor is 
larger. At d = 1.3, it is - 8% for (3,2), (4,2), and (4,4) 
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and ~ 15% for (5,4) and (6,6). This result is qualitatively 
and even quantitatively in good agreement with a result 
in an analytic result presented in Appendix B. 

Calibration for the gravitational wave luminosity 

In Fig. 2, we show the gravitational wave luminosity 
as a function of v 2 . We plot the numerical results (solid 
circles), 2.5PN formula (solid line), 2PN formula (dashed 
line), 1.5PN formula (dotted line), and 1PN formula (dot- 
dashed line). Since v 2 is small in this case, the 2PN and 
2.5PN formulas almost coincide, and the gravitational 
wave luminosity is mostly determined by (2, 2) mode. As 
in the case of the wave amplitude, numerical results agree 
with 2PN and 2.5PN formulas within a small underesti- 
mation by ~ 1.5% for distant orbits. As explained above, 
this error is of numerical origin. For close orbits, the tidal 
effects slightly increase the magnitude beyond the post 
Newtonian formulas, but the amplification is not very 
large (by ~ 5% at d = 1.3). 

Although the effect of the tidal deformation is signif- 
icant for higher multipole components of gravitational 
waves, their contribution to the total luminosity and 
wave amplitude is very small, because the magnitude of 
the (2,2) mode is much larger than others. The amplifi- 
cation factor in the gravitational wave amplitude and lu- 
minosity due to tidal deformation is expected to depend 
strongly on d but weakly on the compactness. Thus, even 
for binaries of large compactness, we expect that the am- 
plification is ~ 2% for the amplitude and ~ 5% for the 
luminosity at the innermost binary orbit, d ~ 1.3. 



3. Effect of location of outer boundary in extracting 
gravitational waves 

As a final calibration, we investigate the effect of outer 
boundary conditions on gravitational wave amplitudes, 
because the outer boundaries are imposed at a finite ra- 
dius. In Fig. 3, we plot the wave amplitude for the 
(2,2) mode as a function of r/A in the case d = 1.3 
(v 2 ~ 0.035). We plot two curves. One (solid line) is 
1-^22 M 1/1-^22 (^ = »" max )| which is obtained by imposing 
the outer boundary condition at r — r max = 55A. The 
other is the result for the following experiment; we im- 
pose the outer boundary condition for a wide range of 
the radius as 0.1A < r max < 55A and compute \H22{r — 

rmax) |- In this case, we plot |F 2 a(r = r max )\/H 2 2(r = 
r mal = 55A)|. We find that (1) if we impose the outer 
boundary condition at r > 5A (10A), the wave amplitude 
can be computed within 0.3% (0.1%) error, (2) if we want 
to compute the wave amplitude within 5% error, it is nec- 
essary to choose the outer radius as r max > 1.5A, and (3) 
even if we impose the boundary condition at r max ~ 0.6A, 
the wave amplitude can be estimated within 15% error. 
In the computation of this paper, we always impose the 
boundary condition at r > 15A, implying that the numer- 
ical error of the wave amplitude associated with the loca- 



tion of the outer boundaries is negligible (much smaller 
than other numerical errors). 

An interesting finding is that even if we imposed the 
boundary condition in the local wave zone (or in the dis- 
tant near zone) at r max ~ A, the wave amplitude could 
be estimated only with a ~ 10% error. In our recent 
simulation on the merger of binary neutron stars, the 
outer boundaries are located in a distant near zone or 
in a local wave zone \r ~ (0.6 — 2)A depending on the 
stage of the merger] pq] . The present results indicate 
that even with this approximate treatment of the outer 
boundary conditions, the gravitational wave amplitude 
could be computed within about a 10% error. 



C. Gravitational waves from compact binaries 

Next, we perform a numerical computation, adopting 
more compact neutron stars. According to models of 
spherical neutron stars, the circumferential radius of re- 
alistic neutron stars of mass M g = 1.4M Q where M Q de- 
notes the solar mass is in the range between ~ 10km and 
~ 15km. This implies that the compactness (M/R)^ 
is in the range between ~ 0.14 and ~ 0.21. Thus, we 
choose (M/R)oo — 0.14 and 0.19 as examples (cf. Table 
I for M g and Mo and Table III for the relevant quantities 
of the quasiequilibrium sequences). 

In Figs. 4-7, we plot the total energy E t and the angu- 
lar momentum J as a function of v 2 for (M/R)^ = 0.14 
and 0.19. They are normalized by Mo and 4Mq to be 
nondimensional. For comparison, we also plot the en- 
ergy and angular momentum for binaries of nonspinning 
stars derived in the 2PN approximation fj as 



E 2 pn = -r)M g v 2 1 - -—r+v 2 - 



^2PN 



9 + r) 2 81 - 57r/ + if 

+2E b , 

r]M g v 2 ( 9 + ?7 2 2(81 — 5777 + ?7 2 ) 



(5.12) 



n 



(5.13) 



where Ef, has to be added in the energy in comparison 
because in E t not only the binding energy between two 
stars but also the binding energy of individual stars is 
included. In 0, J 2 pn is not shown but it is easily com- 
puted from the relation dE = QdJ for the point mass 
case. Figures 4 and 5 show that for distant orbits and 
for (M/i?)oo = 0.14, the numerical results are fitted well 
with 2PN formulas except for a possible small systematic, 
numerical error. This indicates that for mildly relativis- 
tic orbits, higher post Newtonian terms as well as for 
quasiequilibrium binary solutions which we do not take 
into account in this paper are not very important. For 
close orbits as d <■ 1.6, the deviation of numerical results 
from the 2PN formula becomes noticeable. This devi- 
ation seems to be due to the tidal effects because the 
deviation increases rather quickly with increasing v 2 . (If 
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post Newtonian corrections are relevant, the deviation 
should be proportional to a low power of v 2 . On the 
other hand, if tidal effects are relevant, the deviation is 
proportional to d~ 6 oc v 12 p6f.) For (M/R)^ = 0.19 and 
v 2 > 0.1, the coincidence between numerical and 2PN re- 
sults becomes worse even for distant orbits, in particular 
for J. This indicates that effects of third and higher post 
Newtonian corrections could not be negligible for such 
compact binaries. Also, the effects of hij for solutions 
of quasiequilibrium binary neutron stars might not be 
negligible. 

In Figs. 8-11, we show the wave amplitude for the 
(2,2) mode, H22, and the gravitational wave luminosity 
as a function of v 2 for (M/R)^ — 0.14 and 0.19. The am- 
plitude and luminosity are normalized by the quadrupole 
formulas M g v 2 and (2/5)v 10 , respectively. For compari- 
son, we show the 1PN, 1.5PN, 2PN, and 2.5PN formulas. 

v 2 in these sequences of compact binaries is in the 
range between 0.05 and 0.155. The frequency of grav- 
itational waves can be written as 



GW 
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2. 8 M t 
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Thus, if we assume that the total mass of the binary is 
2.8M Q , /gw for binaries presented here is in the range 
between 250Hz and 1350Hz. 

Since convergence of the post Newtonian expansion 
is very slow for v 2 > 0.05, no post Newtonian formu- 
las fit well with numerical results for the whole range of 
v 2 from 0.05 to 0.15. For distant orbits, the numerical 
results agree relatively better with the 2.5PN formulas 
than with lower post Newtonian formulas both for the 
(2,2) mode wave amplitude and for the luminosity. For 
close orbits, on the other hand, the numerical results de- 
viate highly from 2.5PN formulas as well as from other 
formulas. This deviation is due either to the tidal ef- 
fect or to the higher post Newtonian corrections. As we 
show in the small compactness case, the tidal effect could 
amplify the gravitational wave amplitude and luminosity 
by several percent. Therefore, it certainly contributes to 
this deviation. However, the difference between numeri- 
cal results and 2.5PN formulas for v 2 > 0.1 is too large 
to be explained only by the tidal effect. Thus, we con- 
clude that higher post Newtonian corrections affect this 
difference significantly. To explain the behavior of nu- 
merical curves, third or higher post Newtonian formulas 
are obviously necessary . The magnitude of the error 
associated with the neglect of will be estimated in 
Sec. V E. 



D. Validity of assumption for quasiequilibrium 

In this paper, we have assumed that the orbits are in 
quasiequilibrium. As we define in Sec. I, the assumption 
is valid only in the case when the coalescence timescale is 
longer than the orbital period. Here, we assess whether 



the assumption holds for close orbits. To estimate the 
coalescence timescale, we compute 
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dE t 



{-dE/dt) d{v 2 ) 



d(v 2 ), 
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where vq denotes v at an innermost stable orbit, v 2 
should be taken as v 2 at the innermost stable circular 
orbit (ISCO, i.e., the minima for E t and J as a function 
of v or £1) of binaries. However, for irrotational binary 
neutron stars of identical mass with n = 1, the ISCO does 
not exist. As we discussed in |l7j , two neutron stars could 
start mass transfer from their inner edges for d < 1.25, re- 
sulting possibly in a dumbbell-like structure of two cores. 
Even if the shape varies, however, the energy and angu- 
lar momentum are likely to continuously decrease with 
decreasing separation between two cores for d < 1.25, 
and their quasiequilibrium states are mainly determined 
under the influence of general relativistic gravity and the 
tidal interaction between the two cores. Thus, we use an 
extrapolation for the computation of E t and dE/dt for 
d < 1.25 using data points for d > 1.25. A fitting formula 
for E t is constructed using the data points at d = 1.25, 
1.3, 1.4, 1.5, 1.6, and 1.8 as 



Et — &0 + a l v + a 2 v + a 3 v + a 6"V 



1:2 
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where the last term denotes the effect of a tidal deforma- 
tion ||(| . For the fitting, we use the least squares method. 
In Figs. 12 and 13, we show E t in the fitting formula as 
a function of v 2 for (M/R)^ = 0.14 and 0.19. It is found 
that the energy curves around the innermost binary or- 
bit (at d = 1.25) are well fitted by this method and that 
the minimum of the energy appears. We define v 2 as the 
value at the minimum. This minimum is induced by the 
last term of Eq. (5.16) for the case of moderately large 
compactness as (Af/i?) 00 = 0.14. For large compactness 
as (M/i?)oo = 0.19, the minimum appears even without 
the term associated with the tidal interaction, and with 
the tidal term, v 2 at the minimum becomes smaller than 
that without the tidal term. This indicates that not only 
the tidal term but also general relativistic gravity plays 
a role for determining the minimum for such compact 
binaries. 

From Figs. 10 and 11, the gravitational wave luminos- 
ity near the innermost binary orbit at d ~ 1.25 [v 2 ^0.11 
for (M/i?)oo = 0.14 and v 2 ~ 0.15 for (M/R)^ = 0.19] 
may be approximated by (2/5)Ct> 10 where C is a con- 
stant, ~ 0.85 for (M/-R)oo = 0.14, and - 0.80 for 
(A//i?)oo = 0.19. Hence, we use this simple formula 
for the luminosity instead of detailed extrapolation for 
d< 1.25. 

One may think that this procedure is too rough. 
However, it would be acceptable because the evolution 
timescale from the innermost binary orbit at d — 1.25 to 
the minimum found from the fitting formula is ~ 1/3 and 
~ 1/10 of the orbital period at d = 1.25 for (M/R)^ = 
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0.14 and 0.19, respectively. Thus, this rough estimation 
does not cause any serious numerical error. (In other 
words, the orbit at d = 1.25 is close to the ISCO for both 
compactnesses.) 

In Fig. 14, we show i coa i as a function of v 2 for 
(M/i?)oo = 0.14 (solid circles) and 0.19 (solid squares). 
For comparison, we plot the orbital period (solid lines) 
and coalescence time for the Newtonian binary of two 
point masses (i.e., 5M t /(64v 8 ); see All the quantities 
are plotted in units of M t = 2M g . The coalescence time 
becomes equal to the orbital period at d = o? cr i t ~ 1.4 
and v 2 ~ 0.10 for (M/R)^ = 0.14 and at d = d clit ~ 1.7 
and v 2 — 0.125 for (M/R)^ = 0.19. As mentioned in 
Sec. I, assuming the quasiequilibrium state for binary 
neutron stars is appropriate only for distant orbits as 
d > dcrit- At d ~ d cr it, it is likely that an adiabatic cir- 
cular orbit gradually changes to a plunging noncircular 
orbit. This implies that the quasiequilibrium treatment 
for close binary neutron stars can introduce a certain sys- 
tematic error, although it seems still to be an adequate 
approximation as long as the radial approaching velocity 
is much smaller than the orbital velocity (see below). 

The coalescence time we derived here is much shorter 
than the Newtonian coalescence time of two point masses 
for close orbits, although the two results are in better 
agreement for v — > 0. The main reason for the disagree- 
ment is that the variation of the curve for E t becomes 
very slow due to tidal effects for close orbits. [Recall that 
the coalescence time depends strongly on dE t /d(v 2 ).] In 
the absence of the tidal effect, the shape of the curve for 
E t would be similar to that for a binary of point masses, 
so that variation of the energy near the innermost bi- 
nary orbit would not be very slow, and consequently the 
coalescence time would not be as short as the above nu- 
merical results. 

In Fig. 15, we show the ratio of an average, rela- 
tive radial velocity between two stars [defined as i£ ve = 
2Rod(d) /dt] to an orbital velocity v as 



2i? d(d) _ 1 / dE\/ dE t Y 
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The solid and dashed lines denote the numerical results 
for (M/R)oo = 0.19 and 0.14. The dotted line denotes 
the Newtonian result for two point masses (i.e., 16v 5 /5; 
see ||). Figure 15 shows that at d = d CT it, the radial 
velocity is still ~ 2% of the orbital velocity, but it be- 
comes ~ 10% of the orbital velocity near d = 1.25. It 
is also found that the Newtonian formula underestimates 
the radial velocity by several 10% for orbits at d = d C rit- 
For (M/i2)oo = 0.14, the factor of this underestimation 
is rather large, because in this case, the tidal effect which 
increases the radial velocity is significant at d = d CI n. 

It is appropriate to give the following word of caution. 
Since assuming quasiequilibrium states for binary neu- 
tron stars is not very good for d < d cr it, the velocity 
ratio derived for such close orbits might not be a good 



indicator. For d < d cr itj the orbits of a binary could 
deviate from the equilibrium sequence derived in this 
paper. The radial velocity computed in this paper de- 
pends strongly on [dE t /d(d)]^ 1 which becomes very large 
around d ~ 1.25. For a real evolution of binary neutron 
stars, the time evolution of E t could be fairly different 
from the curve for the quasiequilibrium sequence. To 
derive the radial velocity appropriately, numerical sim- 
ulation with an initial condition at d ~ <i cr it may be a 
unique method for this final phase. 



E. hij in the near zone 

In this paper, we have computed quasiequilibrium 
states assuming that the three-metric is conformally flat. 
For the computation of gravitational waves, we also 
adopt a linear approximation in hij, assuming that the 
magnitude of hij is much smaller than unity. In this sec- 
tion, we investigate whether these assumptions are in- 
deed acceptable even for close and compact binaries of 
neutron stars. In the following, we compute the near- 
zone metric of (2,0) and (2,2) modes because they are 
the dominant terms. 

In Figs. 16 and 17, we show h rr and h w computed 
from (2,0) and (2,2) modes along the axis which connects 
the mass centers of two stars for (M/R)^ = 0.14 and 0.19 
and for d = 1.3. For comparison, we also show -00 — 1. 
The centers of the two stars are located at r ~ 0.05A. 
It is found that the magnitude of each mode is < 0.1 
and sufficiently smaller than — 1, which denotes the 
deviation from flat space in the conformal part of the 
three-metric. Second post Newtonian studies |3^j3l[] in- 
dicate that hij is a quantity of 0(v A ) and of 0(v 2 ) smaller 
than 4(V>o — 1), and the numerical results here agree ap- 
proximately with the post Newtonian results. Since the 
magnitude of hij is smaller than 0.1 even for strongly 
relativistic cases, neglecting the nonlinear terms of hij 
appears to be acceptable as long as we allow an error of 
^1%. However, the magnitude of hij is not small enough 
to neglect the linear term. Thus, quasiequilibrium states 
computed in the conformal flatness approximation likely 
contain a systematic error of certain magnitude. 

From a simple order estimate using basic equations for 
computation of quasiequilibrium states, several quanti- 
ties could be modified in the presence of hij as 
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where quantities with 5 denote the deviation due to the 
presence of hij . 

For (M/R)oa = 0.19 and for close orbits as d = 1.3, 
the absolute magnitude of hij at the location of stars is 
typically ~ 0.05. This implies that neglecting hij might 
induce a systematic error of O(10~ 2 ) for and J and of 
0(1O -3 ) for p, ijj, and M for close and compact binaries. 
These systematic errors might also induce a systematic 
error for the frequency and amplitude of the gravitational 
radiation of O(10~ 2 ). Obviously, hij cannot be neglected 
for close and compact binaries if we require an accuracy 
within a 1% error. 

For (M/Ffyoo — 0.14 and d — 1.3, the magnitude of hij 
is about half of that for (M/R)^ = 0.19, i.e., - 0.02 at 
the location of stars. This is reasonable because hij is of 
0(i> 4 ). Thus, for smaller (M/R)^, the conformal flatness 
approximation becomes more acceptable. However, even 
for (M/iJjoo = 0.14, the magnitude of the systematic 
error due to the neglect of h^ could be larger than 1% for 
close orbits, implying that it seems to be still necessary 
even for neutron stars of mildly large compactness to take 
into account hij to guarantee an accuracy within a 1% 
error. 

Finally, we carry out an experiment: In solving equa- 
tions for the nonaxisymmetric part of hij, we have im- 
posed an outgoing wave boundary condition since it 
obeys wave equations. This boundary condition is nec- 
essary to compute gravitational waves in the wave zone. 
However, to compute the near-zone metric for r <C A, the 
term (£ k dk) 2 hij in the wave equation is not very impor- 
tant because its magnitude ~ hij/X 2 is much smaller 
than that of A/iy ~ hij/d 2 or hij/R 2 . This implies 
that solving an elliptic type-equation, neglecting the term 
(£ h dk) 2 hij, could yield an approximate solution for hij 
in the near zone. In this experiment, thus, we solve the 
elliptic-type equation for A22 as an example and compare 
the results with that obtained from the wave equation to 
demonstrate that this method is indeed acceptable for 
computation of the near-zone metric. 

The elliptic-type equation for A22 is solved under the 
boundary conditions 



at r = 0, and 



dA x 
dr 



.4, 



= 0, 



(5.23) 
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at r ^S> A. The outer boundary condition is determined 
from the asymptotic behavior of the source term. 

In Fig. 18, we show h rr and h vv /r 2 computed from 
two different equations of different asymptotic behaviors 
in the case when (M/R)^ = 0.19, d = 1.3, and v 2 ~ 0.15 
(i.e., in the highly relativistic case). Note that the centers 
of stars are located at r ~ 0.052A and the stellar radius is 
0.040A. It is found that the two results agree fairly well 
for r < 0.1 A where the stars are located. The typical 



magnitude of the difference between the two results is of 
O(10" 3 ). 

Accor ding to a post Newtonian theory in the 3+1 for- 
malism |37|[$1|, the difference between the two results 
denotes a radiation reaction potential of 2.5PN order. In 
our present gauge condition, the 2.5PN radiation reac- 



tion potential is written as 37 

h R = 4rf 3 ^, 
y 5 dt 3 



(5.25) 



where -tij denotes the trace-free quadrupole moment. 
For Newtonian binaries of two point masses in circular 
orbits in the equatorial plane, we find 



h t = ~ h yy = ~i v ' sin *' 



h* y = -w 5 cos*, 



(5.26) 



where other components are vanishing. Equation (5.2£) 
indicates that the magnitude of hfj is of O(10 -3 ) even 
for v 2 ~ 0.1. Therefore, the results shown in Fig. 18 are 
consistent with the post Newtonian analysis. 

As mentioned above, the configuration of binary neu- 
tron stars and the orbital velocity are determined by 
quantities in the near zone. Thus, for obtaining a re- 
alistic binary configuration and orbital velocity taking 
into account hij , solving modified elliptic- type equations 
instead of the wave equations for hij may be a promising 
approach. 



VI. SUMMARY AND DISCUSSION 

We present an approximate method for the compu- 
tation of gravitational waves from close binary neutron 
stars in quasiequilibrium circular orbits. In this method, 
we divide the procedure into two steps. In the first 
step, we compute binary neutron stars in quasiequilib- 
rium circular orbits, adopting a modified formalism for 
the Einstein equation in which gravitational waves are 
neglected. In the next step, gravitational waves are com- 
puted solving linear equations for hij in the background 
spacetimes of quasiequilibria obtained in the first step. 
In this framework, gravitational waves are computed by 
simply solving ODEs. The numerical analysis in this pa- 
per demonstrates that this method can yield an accurate 
approximate solution for the waveforms and luminosity 
of gravitational waves even for close orbits just before 
merger in which the tidal deformation and general rela- 
tivistic effects are likely to be important. 

From numerical results, we find that tidal and gen- 
eral relativistic effects are important for gravitational 
waves from close binary neutron stars with d < 1.5 and 
v 2 > 0.1. As a result of tidal deformation effects, the 
amplitude and luminosity of gravitational waves seem to 
be increased by a factor of several percent. It is also indi- 
cated that convergence of the post Newtonian expansion 
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is so slow that even the 2.5PN formula for the luminosity 
of gravitational waves is not accurate enough for close 
binary neutron stars of v 2 > 0.1. 

In Sec. V E., we indicate that the magnitude of a sys- 
tematic error in quasiequilibrium states associated with 
the conformal flatness approximation with hij = is 
fairly large for close and compact binary neutron stars. 
To investigate the quasiequilibrium states and associ- 
ated gravitational waves more accurately, we obviously 
need to improve the formulation for gravitational fields 
of quasiequilibrium states. Thus, in the rest of this sec- 
tion, we discuss possible new formulations in which an 
accurate computation will be feasible. Although a few 
strategies have been already proposed Jl(],[ll| , there seem 
to be many other possibilities, as we here propose some 
different methods in the case when we assume the pres- 
ence of the helical Killing vector. 

The most rigorous direction is to solve the full set of 
equations derived in Sec. II. However, to adopt this, we 
have to resolve several problems. One of the most serious 
problems is that the total ADM mass diverges because 
of the presence of standing gravitational waves in the 
whole spacetimes. This implies that the spacetime is not 
asymptotically flat, and it appears that we have to im- 
pose certain outer boundary conditions in the local wave 
zone just outside the near zone (i.e., at r ~ A). In this 
case, it is not clear at all what the appropriate boundary 
condition is for geometric variables. As we indicated in 
Sec. IV, if we impose an inappropriate outgoing wave 
boundary condition in the local wave zone, the error in 
the gravitational wave amplitude could be rather large. 
Thus, for adopting this strategy, we need to develop ap- 
propriate outer boundary conditions for the gravitation 
fields. We emphasize that numerical computation with 
rough boundary conditions leads to a fairly inaccurate 
numerical result in this strategy. 

One of strategies for escaping this "standing wave 
problem" is to adopt a linear approximation with respect 
to dthij. Note that the divergence of the ADM mass and 
related problems for imposing outer boundary conditions 
are caused by the terms AijA 13 in the equations for a and 
tp and by the term AikA^ in the equation for hij which 
contain the quadratic terms of dthij and hence behave as 
0(r~ 2 ) in the wave zone. Thus, if we neglect the nonlin- 
ear terms of dthij in the equations of a, tp, and hij, there 
is no problem in solving these equations with asymptot- 
ically flat outer boundary conditions in the distant wave 
zone. As we indicated in Sec. V E, nonlinear terms of 
hij are small in the near zone, so that neglect of them 
would not cause any serious systematic error. The ne- 
glect is significant in the wave zone because it changes 
the spacetime structure drastically. However, as men- 
tioned in Sec. IV, this linearization may be considered 
as a prescription to exclude the unphysical pathology as- 
sociated with the existence of the standing wave. One 
concern in this procedure is that the solutions derived in 
this formalism do not satisfy the Hamiltonian constraint 



equation, because we modify it, neglecting the nonlinear 
terms of dthij. However, as long as the magnitude of the 
violation is smaller than an acceptable numerical error, 
say, ~ 0.1%, this method would be acceptable. 

Even simpler method is to change the wave equation 
for h^ to an elliptic-type equation, neglecting the term 
(£ k dk) 2 hij. By this treatment, we can exclude the prob- 
lem associated with the existence of standing waves. In 
this case, we do not have to neglect nonlinear terms of 
h^ because they do not cause any serious problems in 
the distant zone. As shown in Sec. V E, even if wc 
solve the elliptic-type equation for hij, the solution in 
the near zone likely coincides well with the solution ob- 
tained from the wave equations. This indicates that this 
treatment could yield an accurate approximate solution 
for the near-zone gravitational field and matter configu- 
ration of binary neutron stars. In this case, gravitational 
waves cannot be simultaneously computed. However, as 
we have shown in this paper, we can compute gravita- 
tional waves in a post-processing. 

The method we should choose depends strongly on our 
purpose. If one would want to obtain an "exact" so- 
lution in the presence of the helical Killing vector, we 
should choose the first one, even though it may be an 
unphysical solution. However, if we would want to ob- 
tain a reasonably accurate, physical solution or to obtain 
theoretical templates of reasonable accuracy, say, within 
0.1% error, some approximate methods such as second 
and third ones may be adopted. We think that our pur- 
pose is not to obtain the unphysical, exact solution but 
to obtain a reasonably accurate physical solution which 
can be used as theoretical templates. In using second 
and third methods, we do not need new computational 
techniques or large-scale simulations. Furthermore, com- 
putational costs will be cheap. For these reasons, we 
consider that the second and third methods are promis- 
ing. 
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APPENDIX A: SOME FUNDAMENTAL 
CALCULATIONS 

With the expansion of hjj in terms of tensor harmon- 



ics functions such as Eq. ( 4.1SQ , the components of the 



Laplacian of hij are written as 



Ah r 



E 



A" 



' A 1 



V + 6 



A, 



4A; 



-Bu 



Y h 
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e Hi m Y h 



Ah r6 = 

l,m 

E 



Xi +4 



-o/m H H *l 



lm 



d e Y 



i rn 



r" Xl+4 r + ~ 2Az n 



sin ( 



= E ( H lm 9 e Y lm + H hrT^fj > 

Z,rn ^ ' 

A/i ry = Y,( H Ld v Y lm - Hf m smed e Yi m ), 



= E 
E 



// 2 , A/ 2 



Xl m 

sin 9Wi m 



l,m 

^ = E(-l H LY lm + HLw lm -H;' x 



l,m 



lm sm9 



Ah 



r 2 sin 2 i 



E (—HLYim - H( m W lm + HL§S) ■ 



(Al) 



APPENDIX B: POST NEWTONIAN 
WAVEFORMS OF GRAVITATIONAL WAVES 
FROM A BINARY OF TWO ELLIPSOIDAL 
STARS 

The leading terms (up to first post Newtonian order) 
in the wave zone expansion of the gravitational waveform 
are decomposed in terms of radiative multipoles as [p8| 



Hj — J ^ I', ,i i 



(2) hi + ^NjVl mkl + ^e mn(k Wj l)m N n 



, (Bi) 



where D is a distance from a source to an observer, 
(")/(<) = <PI/dt n , X (kl) = {X kl +X lk )/2, N k = x k /r, 
Cij k is a completely antisymmetric tensor, and 



Pijkl = {S ik - NiN k )(Sji - NjNi) 

~(5ij ■ - NiNjWki - N k N t ). 



(B2) 



hj, lijkt lijki, Jij, and Jijk in Newtonian order are writ- 
ten in the form 



I\j Qij ^ijQkk, 



(B3) 



lijk — Qijk — - ySijQkii + SikQju + SjkQiiij , (B4) 

lijkl — Qijkl y ^&ij Q klnn ~t~ &ikQ jinn ~\~ ^ilQjknn 



~t~ &jkQilnn ~\~ ^jlQiknn ~t~ ^klQijr 



+ ( SijS k i + SikSji + 5uSj k ) (),„,,. 



(B5) 
(B6) 
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Sij(2Skii + Suk) 



+ Sik{2Sju + Suj) + 5jk{2Sm + Sm) 



where 



lij — = J P x%x ^ ' ' ' d 3 x, 
Sij...k = / px l x j ■ ■ ■ e k i n x l v n d 3 x. 



(B7) 

(B8) 
(B9) 



Here, we consider irrotational binary neutron stars of 
equal mass in equilibrium circular orbits with angular 
velocity f2 in Newtonian gravity for the computation of 
the above multipolc moments. For simplicity, we assume 
that the shape of each star is ellipsoidal and, in a rotat- 
ing frame, the stars are located along the x axis which 
coincides with the semimajor axis. Defining the coordi- 
nates in the rotating frame as (X, Y, Z) and denoting the 
separation between centers of two stars as 2d, we have 
the following nonzero components for Qij... : 

Qxx = 2(M N d 2 + Qx), Q YY = 2Q 2 , Q zz = 2Q 3 , 
Qxxxx = 2(M N d 4 + 6Qid 2 + Q u ), 
Qxxyy = 2(Q 2 d 2 + Q 12 ), 
Qxxzz = 2(Q 3 d 2 + Q 13 ), 
Qyyyy = 2Q22, Qyyzz = 

Qzzzz - 2Q 33 , (BIO) 

where Mn is the Newtonian mass of one star, and Q k 
and Qki denote the quadrupole moments and 2 4 -pole mo- 
ments of each star (1, 2, and 3 denote xx, yy, and zz, 
respectively). In the inertial frame, the nonzero compo- 
nents of Qij are written as 

Qxx = C 2 QXX + S 2 QYY, Qyy = S 2 QxX + C 2 QyY, 

Qxy = sc(Qxx — Qyy), Qzz — Qzz, 



Qxxxx = c 4 Qxxxx + s a Qyyyy + €>c 2 s 2 Qxxyy, 
s 4 Qxxxx + c 4 Qyyyy + 6c 2 s 2 Qxxyy, 



tyyyy 



exxyy 



c 2 s 2 {Qxxxx + Qyyyy 
+ (c 4 + s 4 -6c 2 s 2 )Qxxyy, 

Qxxxy — C 3 sQxXXX + CS 3 QyYYY 
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- 3cs{c 2 - s 2 )Qxxyy, 

Qxyyy = CS QxXXX + C sQyYYY 

+ 3cs(c 2 - s 2 )Qxxyy, 
Qxxzz — c 2 Qxxzz + s 2 Qyyzz, 

Qyyzz — S 2 QxXZZ + C 2 QyYZZ, 

Qxyzz — cs(Qxxzz — Qyyzz), 

Qzzzz = QzZZZ, 



(Bll) 



where c = cos(fii) and s = sin(Ot). 

To compute Sij.-k, we need the velocity which is for- 
mally obtained after we solve the hydrostatic equations. 
For simplicity, we here assume the following form in the 
rotating frame : 



V x \ ( q{Y 

V Y } = { Qd{X/\X\) + q 2 {X-d) 

V z \ 



(B12) 



where q\ and q 2 are constants which depend on the or- 
bital separation 2c?. In the case of incompressible fluid, 
this becomes a highly accurate approximate solution |}9| . 
Thus, for a star of a stiff equation of state such as neu- 
tron stars, this assumption would be acceptable. In this 
velocity field, all components for Sij are zero, and we 
have the following nonzero components for Sijk', 



S xxz = V(c 2 Q\ x \xx + s 2 Q\ X \yy) 

- qi(c 2 QxxYY + s 2 Qyyyy) 
+ q2(c 2 Qxxxx + s 2 Qxxyy), 

Sxyz = CS{V{Q\ X \XX - Q\X\YY) 

— qi(QxxYY — Qyyyy) 

+ q2(Qxxxx — Qxxyy)}, 
Sxzx = ~C 2 VQ\x\ZZ + S 2 q\QYYZZ - c 2 q2Qxxzz, 

cs{VQ\ X \zz + qiQYYzz + q2Qxxzz}, 



Sxzy — 



V(s 2 Q 

\x\xx + c 2 Q\x\yy) 



qi(s Qxxyy 



c Qyyyy) 
+ q2(s 2 Qxxxx + c 2 Qxxyy), 
S y zx = -cs{VQ\ X \zz + qiQYYzz + qiQxxzz}, 
Syzy = -s 2 VQ\ X \zz + c 2 qiQ Y Yzz - s 2 q 2 Qxxzz, 
Szzz = VQ\x\zz - qiQvvzz + qiQxxzz, (B13) 

where V = (O — 92 and 

Q\x\xx = 2(Mud 3 + dQi), 

Q l x\YY = 2dQ 2 , Q\x\zz = 2dQ z . (B14) 

Using these multipole moments, we can compute the 
waveforms and luminosity of gravitational waves as 



Dh+ = Mi 



-(l+U 2 )cOs(2*y/22 



1 



+ |(l- U 4 )cos(4*)t; 4 /44 



84 v 
1 
6 



Qu 2 + 1) cos(2*)v 4 /42 



(2u z - 1) cos(2*)w 4 / 32 



(B15) 



Dh x =M N 2usin(2^)v 2 / 2 2 



- -u(l - u 2 ) sm(4*)« 4 / 44 
+ ^-u(7u 2 - 5) sm(2*) V 4 /42 
+ ^(3u 2 -l)usin(2^)v 4 f 32 



dE 
~dt 



d 



+ 



/22 + G) {es^ 2 

5 f2 ^ 1280 2 



3969 



567 



where v = 2Qd, and 



/ 22 = (Qxx - Q YY )/(M N d 2 ), 

Iaa = (Qxxxx + Qyyyy — QQxxYY)/(M^d 4 

f±2 = [Qxxxx — Qyyyy 

- G(Qxxzz - Q Y Yzz)]/(M N d 4 ), 

f 3 2 = [V{Q\ x \xx ~ Q\x\yy - 2Q\X\ZZ) 
+ qi{— Qxxyy + Qyyyy — 2Qyyzz) 
+ q2(Qxxxx — Qxxyy — 2Qxxzz)] 

/(M N nd 4 ). 



(Bi6) 



(B17) 



(B18) 
(B19) 

(B20) 



(B21) 



We note that the subscript of fi m indicates the compo- 
nent in the exp ansio n by tensor spherical harmonic func- 
tions as in Eq. ( 4.18 ). Other modes of nonzero m besides 
the above modes are vanishing due to 7r-rotational sym- 
metry and plane symmetry with respect to the equatorial 
plane of the system. 

fi m indicates the amplification factor of the gravita- 
tional wave amplitude due to tidal deformation. For 
I = m = 2, it can be written as 



h 



I 



M N d 2 ' 



(B22) 



where the second term denotes the correction due to the 
tidal deformation. Following |3q], we write Qk as 



h = —M N a k , 
5 



(B23) 



where aj, is the length of semi axes and n n is a con- 
stant depending on equations of state. For incompress- 
ible fluid, K n — 1, and K n is smaller for softer equations 
of state (in the Newtonian ~ 0.65 for n = 1 

p6|). Thus, the amplitude of quadrupole gravitational 
waves is increased by the tidal deformation by a factor 
0.2k„(o 2 — a 2 )/d 2 . Since d has to be larger than a\ and 
K n < 1, the amplification rate is at most 0.2. (Accord- 
ing to M, it is at most ~ 0.14 because a 2 /a\ > 0.5 for 
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di = d.) We note that a\ — a\ is proportional to d~ 5 so 
that the amplification factor rapidly increases with de- 
creasing orbital separation. 

The amplification factors for higher multipoles are 
found to be 

Qn + Q22 - 6Q12 



/ 44 = 1 + 6 



1 + 6 



Qi~Q 2 

M N d 2 
Qi - Qs 



M N d 4 



fan = 1 



M N d 2 

Q11 - Q22 - 6(Qi2 - Q13) 

M N d 4 
3Qi - Q a - 2Q 3 



3^2 Qi - 91Q2 



+ 9i 



Mud 2 

— Ql2 + <?22 



M N fld 2 



2Q 



[2] L. Blanchet, B. R. Iyer, C. M. Will, and A. G. Wiseman, 

Class. Quant. Grav. 13, 575 (1996). 
[3] e.g. L. Blanchet, in Relativistic Gravitation and Grav- 
itational Radiation, eds. J.-P. Lasota and J. -A. Marck 
(Cambridge University Press, 1997), 33: Prog. Theor. 
Phys. 136, 146 (1999). 
(B24) [4] C. M. Will, in proceedings eighth Nishinomiya-Yukawa 
Memorial Symposium on Relativistic Cosmology, edited 
by M. Sasaki (University Academy Press, inc., Tokyo, 
1994), 83; Prog. Theor. Phys. 136, 158 (1999). 
[5] H. Tagoshi, A. Ohashi, and B. Owen, Phys. Rev. D 63, 
044006 (2001). 

[6] Y. Itoh, T. Futamase and H. Asada, Phys. Rev. D 63, 

064038 (2001): T. Futamase, private communication. 
[71 L. Blanchet , G. Faye, B. R. Iyer and B. Joguet, gr- 



(B25) 



q_2- 



M N fld 4 
?n — Q12 — 2Q13 



M N fld 4 



(B26) 



Thus, it is obviously found that the magnitude of the 
tidal effect for /I44 and fa is about 6 times larger than 
that for /22- (Q3 is slightly larger than but roughly 
equal to Q2 for binary of incompressible fluid J3(J.) "6 
times" implies that the amplitude of gravitational waves 
for these multipoles can be several 10% larger than 
that without the tidal deformation. For a rough esti- 
mation of fa, we use the relations for incompressible 
fluid. In this case, both gi/f2 and q 2 /Q are written as 



(a? - a 2 )/(a 2 
becomes 

/32 = 1 



36f. Thus, the amplification factor 



1 



M N d 2 



3Qi -Q 2 - 2Q 3 



-(3Qi - Q 2 



Q1-Q2 
Q1 + Q2 



0{d~% (B27) 



indicating that the magnitude of the tidal effect on fa 
could be about 4-5 times as large as that of / 2 2- 

All these results demonstrate that the effect of tidal 
deformation on the gravitational wave amplitude is more 
important for higher multipole gravitational waves and 
qualitatively agree with the numerical results in Sec. V. 

In more higher multipole modes such as the I = m = 6 
mode, a term such as Qxxxxxx will contribute. It is 
evaluated as M^d 6 + l5Qxxd 4 + 0(d 2 ), and the amplifi- 
cation factor due to the tidal deformation will be about 
15 times larger than that for the I = m = 2 mode. Thus, 
the effect of tidal deformation for close binary neutron 
stars will be even more significant. 
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(M/R) x 


M 


Mg 


0.050 


0.059613 


0.058124 


0.140 


0.14614 


0.13623 


0.190 


0.17506 


0.16000 



TABLE I. Compactness, baryon rest mass, and gravita- 
tional mass for spherical stars in isolation for T = 2. Note that 
for a maximum mass star, (M/R)^ ~ 0.214, Mo — 0.180, and 
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a 


IS 2 




T 
J 


V Hi j lVl 


1 Q 

l.o 


O.Oz4oX 1U 


l.lO / OX 1U 


i.yoUox iu 


-l.DOol X IU 


1 A 

1.4 


o.oozoX 1U 


i i c7o win -1 
1.10 / oX 1U 


i.yoyux iu 


-l.oUo ( X IU 


1.0 


O.U / O / X 1U 


i i coo win - ^ 
L.lDoZ X 1U 


O HQ /I 7vin~^ 

z.Uo4 / X IU 


i ooQfi win - '' 


l.o 


Z. /0©4X 1U 


l.looo X 1U 


z.lzlZX IU 


-l.o41 / X IU 


2.0 


2.5406x10^ 


1.1589xl0 _i 


2.2108xl0" :j 


-1.3912x10"^ 


2.2 


2.3288x10^ 


1.1592xl0 _i 


2.3012x10"^ 


-1.3769x10"^ 


2.4 


2.1441x10^ 


1.1594xl0 _i 


2.3861 x 10"^ 


-1.3810x10"^ 


2.6 


1. 9890x10^ 


1.1596xl0 _i 


2.4754x10"^ 


-1.4485x10"^ 


2.8 


1.8488x10^ 


1.1598xl0 _i 


2.5533xl0" :j 


-1.1064x10"^ 


3.0 


1.7345x10^ 


1.1600xl0 _i 


2.6449x10"^ 


-1.2734x10"^ 



TABLE II. A sequence of irrotational binary neutron stars in quasiequilibrium circular orbits 
with small compactness (M/R)^ — 0.05. 



d 


2 

V 


M 


J 


E t /Mo 


VE/M 


1.25 


l. 073x10 


2.693x10 


6.831x10 


-1.573x10 


-4.732x10 


1.3 


1.060xl0 _i 


2.693x 10~ x 


6.850xl0~ 2 


-1.571xl0 _i 


-4.731 x 10 _b 


1.4 


1 023x10 _i 


2 694x 1 0~ x 


6 910x10~ 2 


-1 566x1 _i 


-4 161 x10 _b 


1.5 


9.798x10^ 


2.695xl0 _1 


6.994xl0~ 2 


-1.559xl0 _i 


-4.982xl0 _b 


1.6 


9.353x10-* 


2.696X10" 1 


7.092x10^ 


-1. 552x10^ 


-4.833xl0" b 


1.8 


8.510x10^ 


2.698xl0 _1 


7.314xl0~ 2 


-1.537xl0 _i 


-5.756xl0" b 


2.0 


7.768x10^ 


2. 700x10^ 


7.553x10^ 


-1. 524x10^ 


-6.715xl0" b 


2.2 


7.129xl0 _: " 


2.702X10" 1 


7.798 xl0~ 2 


-1.511xl0 _i 


-7.482xl0" b 


2.4 


6.573x10^ 


2. 704x10^ 


8.034xl0~ 2 


-1.500xl0 _i 


-7.465xl0 _b 


2.6 


6.101x10^ 


2.705xl0 _1 


8.283xl0~ 2 


-1.491xl0 _i 


-8.097xl0~ b 


2.8 


5.679x10^ 


2.706xl0~ x 


8.506x10^ 


-1.482xl0 _i 


-7.376xl0 _b 


3.0 


5.327x10^ 


2. 707x10^ 


8.763x10^ 


-1. 475x10^ 


-7.453xl0 _b 














1.25 


1.508xl0 _i 


3.151xl0 _1 


8.555x10^ 


-1.998xl0 _i 


-4.702xl0 _b 


1.3 


1.495xl0 _i 


3.152xl0 _1 


8.567x10^ 


-1.997xl0 _i 


-4.143xl0 _b 


1.4 


1.453xl0 _i 


3.152xl0 _1 


8.610x10^ 


-1.993xl0 _i 


-3.315xl0~ b 


1.5 


1.397xl0 _i 


3.153xl0 _1 


8.678x10^ 


-1.987xl0 _i 


-4.445 xl0 _b 


1.6 


1.337xl0 _i 


3.155X10" 1 


8.765x10^ 


-1.979xl0 _i 


-4.121xl0 _b 


1.8 


1. 220x10^ 


3.158xl0 _i 


8.972xl0~ 2 


-1.962xl0 _1 


-5.387xl0" b 


2.0 


1.115xl0 _i 


3.161xl0 _1 


9.205xl0~ 2 


-1.946xl0 _i 


-7.191xl0 _b 


2.2 


1. 024x10^ 


3.163X10" 1 


9.452xl0~ 2 


-1.931xl0 _i 


-8.056xl0~ b 


2.4 


9.454xl0~ 2 


3.166X10" 1 


9.695xl0~ 2 


-1.917xl0 _1 


-8.982xl0" b 


2.6 


8.778x10^ 


3.168X10" 1 


9.953x10^ 


-1. 904x10^ 


-9.684xl0 _b 


2.8 


8.178x10^ 


3.170X10" 1 


1.019xl0 _1 


-1.893xl0 _i 


-7.780xl0 _b 


3.0 


7.668xl0~ 2 


3.172xl0 _i 


1.046xl0 _i 


-1.883xl0 _1 


-8.335xl0^ b 



TABLE III. The same as Table II but for (M/R)^ = 0.14 (upper) and 0.19 (lower). 
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FIG. 1. Relative difference between the numerical results 
and post Newtonian analytic results as a function of v 2 
for several modes of the gravitational wave amplitude for 
(M/R)^ = 0.05. 
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FIG. 4. The total binding energy Et in units of Mo as a 
function of v 2 for (M/R)^ = 0.14 (solid circles). For compar- 
ison, we plot a curve derived from the 2PN formula (dashed 
line). 
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FIG. 2. The gravitational wave luminosity normalized 
by the quadrupole formula 0.4i> 10 as a function of v 2 for 
(M/R)oa = 0.05. The numerical results (solid circles), and 
1PN (dot-dashed line), 1.5PN (dotted line), 2PN (dashed 
line), and 2.5PN (solid line) formulas are shown. The 1.5PN 
formula is very close to the 2PN formula. 
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FIG. 3. |-ff22(r)/-H22(r m ax)] as a function of r for for the 
case we impose boundary conditions at r max = 55A (solid 
line) and | .H22 (rmax)/H 22{r mBX = 55A)| in an experiment 
of varying r max from 0.1A to r max = 55A (dashed line) for 
{M/R)oo = 0.05 and d = 1.3. 
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FIG. 5. The same as Fig. ^ but for the total angular mo- 
mentum divided by (2Mq) 2 . 




FIG. 6. The same as Fig. |but for {M/R)^ = 0.19. 
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FIG. 7. The same as Fig. |but for {M/R)^ = 0.19. 



FIG. 10. The same as Fig. but for (M/R)oo = 0.14. 
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FIG. 8. Amplitude of gravitational waves for the (2,2) 
mode (H22) normalized by M g v 2 as a function of v 2 for 
(M/R)oo — 0.14 (solid circles). For comparison, we plot the 
results for the 1PN (dot-dashed line), 1.5PN (dotted line), 
2PN (dashed line), and 2.5PN (solid line) formulas. 
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FIG. 9. The same as Fig. U but for (M/R)^ = 0.19. 
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FIG. 11. The same Fig. [fi] but for (M/R)^ = 0.19. Even 
for this highly compact case, the 1.5PN formula (dotted line) 
is very close to the 2PN (dashed line) formula. 
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V 2 

FIG. 12. Fitting formula for Et/Mo around the innermost 
binary orbit (at d = 1.25) for (M/R)^, = 0.14 (solid line). 
The 2PN formula (dashed line) and numerical data points 
are also plotted. 



25 




0.12 0.13 0.14 0.15 0.16 0.17 

FIG. 13. The same as Fig. |l| but for (M/R) x = 0.19. 
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FIG. 14. The coalescence time f coa i as a function of v 2 for 
(M/R)^ = 0.14 (solid circles) and for (M/R)^ = 0.19 (solid 
squares). In comparison, the orbital period and coalescence 
time in the Newtonian point mass case (t coa i — 5M t /(64n 8 )) 
are plotted by the solid and thin dotted lines. The time is 
shown in units of M t (— 2M S ). 
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FIG. 15. The ratio of an aver age, r elative radial velocity 
to an orbital velocity [see Eq. ( 5.17 )1 as a function of v 2 
for (M/R)^ = 0.14 (dashed line) and 0.19 (solid line). The 
thin dotted line denotes the Newtonian formula for two point 

masses (v^ vc = 16i> 6 /5). 
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FIG. 16. Behavior of some of metric components in the 
near zone for {M/R) x = 0.14 and d = 1.3. v 2 ~ 0.106 in this 
case. 
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FIG. 17. The same as Fig. |l| but for (M/R)^ = 0.19 and 
d = 1.3. v 2 ~ 0.150 in this case. 




x/X 

FIG. 18. (2,2) modes of hij in the near zone for 
(M/R)oo — 0.19 and d — 1.3. Dotted lines are results ob- 
tained using a nonwave-type outer boundary condition and 
solid lines are results using an outgoing wave boundary con- 
dition. The centers of stars are located at r ~ 0.052A and 
radius of stars is ~ 0.040A in this case. 
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